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Abstract 

Attitude  control  is  a  requirement  for  most  satellites.  Many  schemes  have  been 
devised  over  the  years  including  control  moment  gyros,  reaction  wheels,  spin  stabi¬ 
lization  and  gravity  gradient  stabilization.  For  low  Earth  orbits,  the  Earth’s  atmo¬ 
sphere  can  have  an  affect  on  a  satellite’s  orbit  and  attitude.  Using  the  atmosphere 
to  control  spacecraft  attitude  has  been  researched  in  the  past  however  very  little  re¬ 
search  has  been  done  using  an  active  feedback  control  system  to  maintain  spacecraft 
attitude. 

This  research  effort  examines  the  feasibility  of  using  the  atmosphere  to  actively 
control  a  spacecraft’s  attitude  using  drag  panels.  Several  variables  affect  the  drag 
force,  of  which,  projected  area  is  the  only  variable  that  can  be  changed  easily.  Adding 
controllable  drag  panels  to  a  satellite  gives  the  ability  to  change  the  projected  area  as 
well  as  the  location  of  the  projected  area.  The  result  of  manipulating  the  projected 
areas  is  a  force  that  is  not  aligned  with  the  center  of  gravity,  resulting  in  an  external 
torque  on  the  spacecraft.  Although  these  torques  are  very  small,  on  the  scale  of 
micro-Newton  meters  and  smaller,  over  time  these  torques  can  be  used  to  change 
the  spacecraft’s  attitude. 

A  linear  computer  model  was  created  using  a  proportional  controller.  This 
model  was  used  to  evaluate  the  effectiveness  of  using  drag  panels  for  attitude  control. 
Results  from  the  simulation  show  that  the  spacecraft  can  recover  from  disturbance 
torques  that  may  cause  a  change  in  attitude  very  effectively  especially  at  low  altitudes 
(200-300km).  At  200km,  the  satellite  is  able  to  recover  from  a  disturbance  in  less 
than  one  hour.  As  the  altitude  increases,  these  settling  times  increase  exponentially. 
At  600km  it  takes  approximately  2  weeks  to  stabilize.  Other  factors  affect  the  settling 
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time  such  as  mass,  and  the  geometric  dimensions  of  the  satellite’s  control  panels  and 
control  arms. 
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Satellite  Attitude  Control 


Using  Atmospheric  Drag 


I.  Introduction 

Over  the  last  several  decades  since  Sputnik  became  the  first  artificial  satellite  to 
orbit  the  Earth  in  1957,  satellites  have  become  an  integral  part  of  our  lives.  Satellites 
serve  many  useful  purposes  such  as  communications,  weather  and  remote  sensing  to 
name  a  few.  As  their  missions  differ,  so  do  their  orbits.  Many  of  these  missions 
require  the  spacecraft  to  be  in  a  low  Earth  orbit  (LEO) 

Some  advantages  for  putting  a  spacecraft  in  LEO  are: 

•  Higher  resolution  for  Earth  sensing  satellites. 

•  Allows  for  smaller/lower  mass  payloads. 

•  Less  costly  to  get  to  LEO. 

•  Shorter  orbital  periods  for  rapid  revisits. 

•  No  fuel  needed  for  placement  into  a  graveyard  orbit. 

There  are  also  disadvantages  for  putting  spacecraft  in  LEO: 

•  Atmospheric  drag. 

•  Finite  orbit  lifetime. 

•  Attitude  control  issues. 

•  Additional  fuel  for  orbit /attitude  maintenance. 

The  atmosphere  is  one  of  the  dominating  contributors  to  orbit  and  attitude 
perturbations  for  spacecraft  in  LEO  and  accounts  for  most  of  the  disadvantages. 
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The  atmosphere  is  the  dominating  contributor  to  these  perturbations  from  altitudes 
ranging  from  150-600  km  [9,  page  81].  The  research  presented  here  only  focuses  on 
these  altitudes  where  the  Earth’s  atmosphere  has  a  significant  effect  on  spacecraft. 

Measurements  show  that  atmospheric  density  tends  to  decrease  exponentially 
with  increasing  altitude  and  can  be  easily  modeled  (see  Appendix  C).  Since  the 
atmosphere  doesn’t  end  abruptly,  there  are  fewer  air  molecules  at  LEO  altitudes. 
As  the  spacecraft  orbits  the  Earth,  it  hits  these  molecules  and  loses  some  of  its 
orbital  energy.  Although  the  energy  loss  is  extremely  small,  over  time  the  losses  add 
up  causing  the  orbit  to  decay  thus  leading  to  a  finite  orbital  lifetime.  There  is  no 
way  to  stop  orbit  decay  other  than  having  onboard  thrusters  for  orbit  maintenance. 
On  the  other  hand,  collisions  with  atmospheric  particles  can  affect  the  spacecraft’s 
attitude.  If  a  satellite  has  a  center  of  pressure  that  is  not  in  line  with  it’s  center  of 
mass  or  center  of  gravity  (CG)  a  small  torque  will  be  produced.  If  the  spacecraft 
utilizes  an  attitude  control  system  other  than  gravity  gradient,  it  will  most  likely 
use  thrusters  for  attitude  control  or  for  momentum  dumping  of  reaction  wheels.  The 
torques  caused  by  the  atmosphere  will  increase  the  spacecraft’s  fuel  consumption 
thus  reducing  it’s  life.  Gyros  and  reaction  wheels  can  also  be  prone  to  failure  which 
can  render  a  satellite  useless. 

1.1  Research  Objectives 

The  exponential  decay  of  the  density  of  the  atmosphere  as  a  function  of  increas¬ 
ing  altitude  is  always  figured  into  satellite  design  for  LEO  spacecraft.  Most  attitude 
control  systems  used  to  overcome  perturbations  tend  to  be  complex  and  expensive. 
The  mission  specifies  the  altitude  of  the  orbit,  from  which  orbital  velocity  and  at¬ 
mospheric  density  can  be  calculated.  The  only  variables  we  can  easily  control  are 
the  projected  areas  and  the  location  of  the  projected  area  relative  to  the  spacecraft’s 
center  of  gravity.  By  changing  the  size  and  location  of  the  projected  area,  torques 
can  be  produced  from  the  drag  the  projected  areas  experience.  This  research  looks 
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into  the  feasibility  of  using  drag  panels  to  produce  torques  on  a  spacecraft  to  damp 
oscillations  and  control  the  spacecraft’s  attitude.  Since  the  forces  on  the  drag  panels 
are  very  small,  the  drag  panels  can  be  as  simple  as  having  a  piece  of  foil  mounted  to 
a  stiff  lightweight  frame,  which  will  then  be  connected  to  an  actuator.  An  attitude 
control  system  using  drag  panels  would  be  less  complex  and  less  prone  to  failure. 
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II.  Background 


2.1  Literature  Review 

Atmospheric  drag  has  had  an  effect  on  most  satellites  starting  with  Sputnik. 
Sputnik  1  was  launched  on  4  October  1957  with  an  initial  apogee  of  approximately 
950  km  and  a  perigee  of  approximately  230  km.  At  its  perigee,  the  atmosphere  had 
a  significant  affect  on  the  satellite  which  eventually  circularized  its  orbit  (apogee  was 
approximately  600  km  by  9  December  1957).  Sputnik’s  orbit  eventually  decayed  to 
the  point  of  re-entry  after  92  days  in  orbit  on  4  January  1958.  Since  then,  many 
satellites,  weather  balloons,  and  sounding  rockets  have  been  launched  to  study  the 
atmosphere.  Some  of  the  studies  that  have  been  done  include  both  active  and  passive 
attitude  stabilization  techniques  using  atmospheric  drag.  Most  of  the  research  in 
satellite  aero  stabilization  deals  with  passive  stabilization.  Little  research  has  been 
done  using  active  attitude  stabilization  using  atmospheric  drag. 

2.1.1  Paddlewheel  Satellites.  By  1968,  about  a  dozen  satellites,  such  as 
Explorer  VI  (see  Figure  2.1),  have  been  launched  to  study  solar  wind  and  the  mag¬ 
netosphere.  While  these  satellites  were  in  orbit,  they  revealed  information  about 
the  aerodynamic  interaction  of  air  molecules  with  the  satellite  surfaces.  According 
to  the  data,  air  molecules  experience  nearly  diffuse  reflections.  Maxwell’s  classi¬ 
cal  model  best  approximated  this  interaction.  Maxwell’s  model  accounts  for  both 
specular  and  diffuse  reflections.  Scientific  analyses  of  the  behavior  of  paddlewheel 
satellites  led  to  a  good  approximation  of  the  accommodation  coefficients  from  which 
an  approximated  drag  coefficient  could  be  determined.  [4] 
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Figure  2.1:  Explorer  VI 


2.1.2  Nanosatellite  Passive  Attitude  stabilization.  Psiaki  performed  a 
study  using  passive  attitude  stabilization  for  a  nanosatellite  [5].  This  satellite  was 
a  cubesat  with  dimensions  of  0.1  nr  for  each  side  and  mass  of  1  kg.  This  satellite 
was  designed  to  use  passive  drag  torques  to  stabilize  the  roll,  pitch  and  yaw  axes 
and  provide  magnetic  damping  on  both  the  pitch  and  yaw  axes.  The  satellite  resem¬ 
bles  a  shuttlecock  used  in  badminton  (see  Figure  2.2).  The  conclusion  of  this  study 
determined  that  the  nanosatellite  performed  reasonably  well  in  simulations. 
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sign 

2.1.3  Shuttle  Hitchhiker  Passive  Aero  stabilization.  Another  study  focused 
on  purely  passive  stabilization  with  magnetic  damping.  This  was  a  NASA  study  for 
the  feasibility  for  a  low  cost,  low  weight,  and  long  life  spacecraft  for  the  gravity  and 
magnetic  Earth  surveyor  (GAMES)  mission  .  It  is  very  similar  to  the  satellite  pic¬ 
tured  in  Figure  2.2  as  far  as  size,  mass  and  method  used  for  aero-stabilization.  Both 
use  magnetic  damping  but  the  physical  design  is  different.  The  shuttle  hitchhiker  is 
a  cylindrical  design  vs.  a  cubesat  with  “feathers”  (see  Figure  2.3).  In  this  design, 
the  center  of  pressure  (CP)  is  aft  of  the  CG  which  tends  to  keep  it  pointed  in  the 
direction  tangent  to  the  orbit.  [7,8] 
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Figure  2.3:  (a)  Shuttle  Hitchhiker  schematic(b)  Cross  section  view 


2.1.4  Active  Aerodynamic  Attitude  Stabilization.  Very  little  research  has 
been  done  with  active  stabilization  using  atmospheric  drag.  Ravindran  and  Hughes 
[6]  researched  an  Earth  oriented  satellite  using  drag  panels  for  attitude  control.  Their 
research  was  closely  related  to  the  research  presented  here,  but  their  satellite  design 
was  very  different.  Their  design  consisted  of  long  cylindrical  body  where  the  axis 
through  the  cylinder  is  aligned  with  the  roll  axis  (see  Figure  2.4).  The  panels  are 
turned  so  they  are  in  the  most  streamlined  position  when  no  control  is  necessary. 
For  pitch  and  yaw  controls,  the  appropriate  panels  would  be  rotated  90°  to  increase 
the  drag  which  is  offset  from  the  CG  thus  producing  a  torque.  For  roll  control,  they 
would  all  be  rotated  a  specified  angle  turning  the  spacecraft  into  a  propeller  like 
configuration.  Results  from  their  analysis  are  similar  to  results  obtained  from  this 
research.  [6] 
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Figure  2.4:  Active  Control  Design 


2.2  The  Atmosphere 

The  atmosphere  has  been  studied  for  quite  some  time.  Dozens  of  satellites, 
weather  balloons  and  sounding  rockets  have  been  launched  to  answer  the  question 
“How  does  the  pressure,  temperature  and  density  change  as  altitude  changes?”  An¬ 
swers  to  this  question  have  been  crucial  to  the  study  of  aeronautics  and  astronautics. 
Most  people  believe  that  space  is  empty  and  void.  The  truth  is,  space  is  mostly  void 
and  empty  however  there  are  particles  flying  around.  According  to  the  1976  standard 
atmospheric  model  [3],  the  atmospheric  density  can  be  approximated  by  using  an 
exponential  model.  See  Appendix  C  for  details  on  how  the  density  was  modeled.  The 
models  used  for  atmospheric  density  only  predict  an  average  value  based  on  solar 
activity.  Solar  activity  can  have  a  significant  affect  on  the  actual  density  by  heating 
and  expanding  the  atmosphere  to  higher  altitudes.  Although  the  atmospheric  den¬ 
sity  model  only  approximates  the  density,  the  dynamic  model  should  still  accurately 
predict  the  dynamic  behavior  of  the  linearized  satellite. 
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2.3  Design  Concepts 

As  progress  was  made  with  this  research,  the  satellite  design  was  altered  to 
simplify  the  linearization  process  for  the  torque  equations.  Three  design  iterations 
were  performed  and  the  satellite  evolved  into  the  final  version  as  presented  in  this 
section  as  Design  3.  This  section  summarizes  these  design  variations  and  why  they 
changed. 

2.3.1  Design  1:  4  Drag  Panels  to  control  Roll,  Pitch  and  Yaw.  The  first 
design  had  4  drag  panels,  each  on  a  control  arm  with  two  degrees  of  freedom  for 
each  drag  panel.  This  design  enabled  3  axis  control  about  the  roll  (&i),  pitch  (b2), 
and  yaw  (63)  axes  (See  Figure  2.5).  The  panels  could  be  rotated  into  the  air  stream 
about  the  b2  axis  for  pitch  control  (top  and  bottom  arms)  and  the  63  axis  for  yaw 
control  (left  and  right  arms).  For  roll  control,  all  control  arms  could  be  rotated 
into  the  air  stream  and  each  panel  rotated  about  each  arm  for  left  and  right  roll 
control.  See  Figures  2.6  (a)  and  (b),  and  2.7  (a)  and  (b).  Design  1  is  very  similar 
in  design  as  the  Explorer  VI  but  with  moveable  panels.  One  drawback  is  that  the 
control  actuators  become  more  complex  having  to  rotate  about  two  different  axes. 
This  design  was  abandoned  since  linearizing  the  equations  of  motino  were  became 
unmanageable.  The  nonlinear  equation  derivations  for  this  design  are  included  in 
Appendix  F. 
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Figure  2.5:  Roll  (Si.),  Pitch  (S2)and  Yaw  (S3)  Axes 


(a)  (b) 

Figure  2.6:  (a)  Panels  retracted  for  normal  flight. 

(b)  Panels  extended  and  rotated  for  left  rolling  maneuver. 
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(a)  (b) 


Figure  2.7:  (a)  Top  panel  extended  for  a  pitch  up  maneuver, 

(b)  Left  panel  extended  for  a  left  yaw  maneuver. 
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2.3.2  Design  2:  Ailerons  to  control  Roll,  4  Drag  Panels  to  Control  Pitch  and 
Yaw.  Due  to  the  complexity  of  the  linearization  of  the  torque  equations  in  Design 
1,  a  revision  to  the  original  design  was  made  to  simplify  the  equations.  This  design 
was  similar  to  the  first  design.  Pitch  and  yaw  were  controlled  using  the  same  drag 
panels.  In  this  case,  the  control  panels  only  had  one  degree  of  freedom  and  now 
only  control  pitch  and  yaw.  “Ailerons”  for  lack  of  a  better  term,  were  added  for  roll 
control.  If  an  actual  satellite  were  to  be  built  and  tested,  either  this  design  or  the 
previous  design  would  be  the  best  choice.  Figure  2.8  shows  the  addition  of  ailerons 
to  the  sides  of  the  spacecraft  bus. 


Figure  2.8:  Satellite  shown  with  ailerons  in  a  right  rolling 

maneuver. 

2.3.3  Design  3:  The  third  and  final  design  was  very  similar  to  the  previous 
design  using  ailerons  for  roll  control.  The  only  difference  between  the  two  models 
is  that  the  ailerons  were  changed  from  flat  plates  to  wedge  shapes.  This  simplified 
the  linearization  by  eliminating  negative  angles  for  the  aileron  surfaces.  This  design 
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is  shown  in  Figure  2.9.  All  simulations  presented  in  this  research  are  based  on  the 
wedge  shaped  aileron  design. 


Figure  2.9:  Satellite  shown  with  wedges  for  ailerons. 
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III.  Methodology 

The  equations  of  motion  and  torque  equations  depend  on  a  satellite’s  geometry, 
physical  dimensions  and  orbit.  When  deriving  the  equations,  all  terms  were  kept  as 
variables  to  give  the  model  more  flexibility  while  running  simulations. 

3.1  Coordinate  Systems 

Three  reference  frames  are  used  for  the  rotational  kinematic  equations.  The 
UB  frame”,  a  body  fixed  reference  frame  is  aligned  with  the  principal  moments  of 
inertia,  or  the  principal  axes  for  the  satellite,  where  b±  is  the  roll  axis,  62  is  the  pitch 
axis  and  63  is  the  yaw  axis.  The  “A  frame”  is  the  local  vertical,  local  horizontal 
(LVLH)  frame  aligned  with  the  orbit  where  hi  is  along  the  orbit  direction,  a 2  is 
perpendicular  to  the  orbit  plane  and  a 3  is  toward  the  center  of  the  Earth.  Finally 
the  inertial  “I  frame”  is  an  Earth  fixed  frame  with  its  origin  at  the  center  of  the 
Earth.  (See  Figure  3.1)  [10,  page  366]. 
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Figure  3.1:  Coordinate  frames  for  a  rigid  body  in  circular 

orbit 

3.2  Rotation  Matrices  Between  Coordinate  Frames 

The  first  step  in  deriving  the  equations  of  motion  is  to  derive  the  rotation 
matrices.  A  3-2-1  (yaw,  pitch,  roll)  rotation  sequence  was  used  (see  Figure  3.2) 
where  each  elementary  rotation  has  the  following  rotation  matrices  associated  with 
it.  Note  c(6i)=  cos(0j),  s(0j)=  sin(0j),  and  6i  are  the  rotation  angles  for  each  of  the 
roll,  pitch  and  yaw  axes  [10,  page  309]. 
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Figure  3.2:  Euler  Angle  Rotations. 


C3(93) 


c(03)  s(03)  0 

-s(93)  c(93)  0 
0  0  1 


C2(92 ) 


c(02)  0  -s(02) 
0  1  0 
s(92)  0  c(02) 


Ci^) 


1  0  0 
0  c(9 i)  s(0i) 

0  -5(01)  c(9\) 


(3.1a) 


(3.1b) 


(3.1c) 
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The  rotation  matrix  to  the  B  frame  from  the  A  frame  is  defined  as  RBA  =  C\{6\)C2(62)C3{63) 
and  when  expanded  becomes  [10,  page  311]: 


Rba 


c{92)c{93)  c(92)s(93)  -s{92) 

s(0i)s(02)c(03)  -  <o i)s(03)  s{e3)s{e2)s{e3)  +  c^yy 3)  s(91y(92) 

c(6i)s(02)c(63)  +  s(^i)s(03)  c{6i)s{02)s{93)  —  s(91)c(93)  c(9i)c(92) 


(3.2) 

To  get  the  rotation  matrix  to  the  A  frame  from  the  B  frame  the  transpose  must  be 
taken  since  RAB=  ( RBA)~ 1  =  (RBA)T ,  [10,  page  311].  RAB  becomes: 


Rab 


c(e2)c(e3)  s(0i)s(02)c(03)  -  cieMh)  c{e1  )s(e2y(o3)  +  s(91y(93) 
c(62)s(63)  s(6i)s(62)s(63)  +  c(9\)c(93)  c(9i)s(92)s(93)  —  s(9i)c(03 ) 
—s(02)  s(9i)c(92)  c(9i)c(92) 

(3.3) 


Finally,  to  get  to  the  B  frame  from  the  A  frame,  or  to  the  A  frame  from  the  B  frame, 
the  equations  are  as  follows  [10,  pages  365-366]. 


k 

a  i 

b2 

=  Rba 

a2 

(3.4) 

_k_ 

a3 

CL\ 

k 

a2 

=  Rab 

b2 

(3.5) 

a3 

b3 

Since  a*  contains  only  orthonormal  unit  vectors,  a  ■  a'  becomes  the  identity 
matrix  and  equation  3.4  can  be  solved  for  RBA: 


k 

biOLi  bia2  b3a3 

rba  = 

b2 

S. 

CL\  CI2 

— 

b2a\  b2a2  b2a3 

b3a\  b3a2  b3a3 

(3.6) 
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3.3  Equations  of  Motion  of  a  Rigid  Body  in  Circular  Orbit 

The  next  step  was  to  derive  the  equations  of  motion  (EOM).  This  derivation  is 
extracted  from  Wie  [10,  pages  365-369].  Several  assumptions  were  made  to  simplify 
the  problem  and  these  assumptions  would  have  little  impact  on  the  problem. 


•  Spacecraft  is  in  a  circular  equatorial  orbit  (constant  orbital  rate). 

•  Uniform  gravitational  field,  (gravity  gradient  torques  are  neglected). 

•  Spacecraft  is  in  low  Earth  orbit  (200-600  km) 

•  Atmospheric  density  is  taken  as  an  average  over  the  orbit  (no  fluctuations  in 
density) . 


3.3.1  Angular  Velocity  Vector.  From  the  coordinate  system  shown  in  Fig¬ 
ure  3.1,  the  angular  velocity  between  the  B  frame  and  the  Earth  frame  /  frame) 
is: 

cu  =  0B1  =  ujba  +  (dAI  (3.7) 

and  written  in  the  body  frame 


id  —  Cdibi  +  id2^2  +  id ^  S3  — 


bi  bo  by 


uq 

id2 

ids 


(3.8) 


Where  uBI  is  the  angular  velocity  vector  of  the  Earth  fixed  frame  to  the  body 
frame,  0AI  is  the  angular  velocity  of  the  I  frame  to  the  A  frame,  and  <2BA  is  the 
angular  velocity  from  the  A  frame  to  the  B  frame. 

Since  the  spacecraft  is  assumed  to  be  in  a  circular  orbit,  uA1  =  — naq,  and 
therefore  equation  3.7  becomes: 


id 


BI 


=  udBA  —  nd2 


(3.9) 
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where  n  is  the  constant  orbital  rate  due  to  a  circular  orbit. 


From  equation  3.5,  a2  can  be  written  as: 


a,2  — 


bi 


c(e2)s(d3) 

s(e1)s(e2)s(e3)  +  c(91)c(03) 

c(91)s(92)s(93)  -  s(91)c(93) 


(3.10) 


To  hnd  ljba,  a  yaw,  pitch,  roll  or  3-2-1  sequence  was  used.  The  first  rotation  about 
the  yaw  (a3)  axis  goes  from  the  A  axis  to  the  A'  axis.  The  second  rotation  about  the 
a'2  axis  takes  us  from  the  A'  to  the  A"  axis  and  finally  the  third  rotation  about  the  a![ 
axis  goes  from  the  A "  to  the  B  frame  (See  Fig.  3.2).  From  the  rotations,  the  angular 
velocity  vector  from  the  A  frame  to  the  B  frame  becomes  [10,  pages  324-326]: 


lua  ‘4  —  93a3  —  93a3 

(3,11) 

uA"A'  =  92a2  =  92a2 

(3.12) 

ujra"  =  9\a[  =  9\b\ 

(3.13) 

Where  6t  is  the  angular  rate  of  the  roll  (i=l),  pitch  (i=2)  and  yaw  (i=3)  axes.  The 
total  angular  velocity  vector  then  becomes: 


ivBA  =  uBA"  +  lJa"a'  +  lja'a  =  93a3  +  92o,2  +  9\b\ 
By  substituting  equations  3.11,  3.12,  and  3.13  into  3.14,  we  get: 


CuBA  = 


bi  b2  b3 


9i 

0 

0 

+  Ci{9i) 

92 

0 

0 

Ci(6,i)c2(6'2) 


0 

0 

9-s 


(3.14) 


(3.15) 
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By  combining  equations  3.7  and  3.8,  and  substituting  in  equations  3.1b,  3.1c,  3.10, 


and  3.15,  the  vector 


bi  b2  b3 


cancels  and  the  equation  becomes  [10,  page  368]: 


uq 

1  0  -s(02) 

0i 

c(02)s(03) 

u2 

= 

0  c(0 i)  s(di)c(02) 

02 

—  n 

s{0l)s(02)s{d3)  +  c(0\)c(0 3) 

(U3 

0  — s(6li)  c(di)c(02) 

03_ 

c(0i)s(02)s(03)  —  s(0i)c(03) 

(3.16) 


By  solving  for  6\,  the  result  is  the  kinematic  differential  equation  for  an  orbiting 
rigid  body  [10,  page  368]: 


0i 

02 

03 


1 

c(02) 

c(02)  s(0\)s(02) 

0  c(0i)c(02) 

c(0l)s(02) 

-s(0i)c(02) 

UJx 

UJ2 

n 

+  <e 2) 

s(03) 

C(02)C(03) 

0  s(0i) 

c(0 1) 

UJ3 

S(02)S(03) 

(3.17) 


The  following  equation  is  the  equation  3.16  expanded: 


Ml 

01  -  s{02)0 3  +  c(02)s(03)n 

U)2 

= 

c{0 1)02  +  s(0\ )c{02)03  +  ns(0i)s(02)s(^3)  +  ^c(6,1)c(03) 

(U3 

-s(0 1)02  +  c(0i)c(02)03  +  nc(0i)s{02) s(03)  -  ns(0l)c(03) 

(3.18) 


Since  0,.  is  very  small,  equation  3.18  can  be  linearized  by  letting  cos (0i)  =  1, 
sin(0j)  =  0i,  sin(#i)  sin(0j)  ~  0  and  sin [6i)6i  ~  0.  The  linearized  equation  becomes 
[10,  page  369]. 


UJl 

1 

1 

CO 

_ 1 

UJ2 

= 

02-n 

u>3 

03  +  n0 1 

(3.19) 
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Next  we  take  the  derivative  of  equation  3.19  to  get  uV 


(hi 

1 

1 

s 

CO 

_ 1 

UJ  2 

= 

02 

UJ3 

O3  +  ndi 

(3.20) 


3.3.2  Time  Derivatives  of  the  Angular  Momentum  Vector.  This  derivation 
follows  Wie  [10,  pages  340-342],  A  rigid  body  has  an  angular  momentum  associated 
with  it,  which  is  expressed  as  H  =  J  ■  ujbi .  By  taking  the  derivative  of  H  the 
rotational  equation  motion  can  be  found  where: 


+  UJ 


BI 


x  H 


(3.21) 


By  using  the  transport  theorem  equation  3.21  becomes: 


J  ■  UJ  +  UJ  X  J  ■  UJ 


M 


(3.22) 


where  J  is  the  moment  of  inertia  matrix  (MOI),  M  is  the  external  moment  vector, 
(jj  is  the  angular  rate  vector  and  uj  is  the  angular  acceleration  vector.  Equation  3.22 
in  matrix  form  becomes: 
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co 
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CN 

"A 

chi 

0  —  (U3  U>2 

J21  J22  J23 
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UJ  3  0  — U>i 
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— UJ 2  UJ\  0 
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J21 
J31 


J 12 
J22 
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Jl  3 
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Mi 
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(3.23) 
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Since  it  is  assumed  that  the  B  frame  is  aligned  to  the  principal  axes,  equation  3.23 
becomes: 


A 

0 

0 

(hi 

0 

—td3 

CJ2 

A 

0 

0 

Wl 

M1 

0 

B 

0 

i02 

+ 

i03 
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= 

m2 
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c 

U>3 

—UJ2 

UJl 

0 

0 

0 

c 

UJ3 

m3_ 

(3.24) 


Where  A  is  the  MOI  about  the  roll  axis,  B  is  the  MOI  about  the  pitch  axis,  C  is  the 
MOI  about  the  yaw  axis,  and  when  multiplied  out  becomes: 


ACj\  —  (5  —  C^)(jJ2UJ3 

M1 

B0J2  —  {C  —  A^u\(jj3 

= 

m2 

C(jj3  —  (A  —  B^(jJ\U2 

m3_ 

(3.25) 


By  substituting  equations  3.19  and  3.20  into  equation  3.25,  linearizing  and  simplify¬ 
ing,  the  equation  becomes: 


ABX  +  (B  —  A  —  C)n03  +  (B  -  C)n261 

Mi 

B02 

= 

m2 

C03  +  (C  +  A  -  B)n0!  +  (B-  A)n263 

m3_ 

(3.26) 


Finally  by  solving  for  6  and  putting  into  matrix  form,  the  equation  becomes: 
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3-4  External  Torques  from  Atmospheric  Drag  (Nonlinear) 


Satellites  in  LEO  experience  an  aerodynamic  drag  force  given  by  [9,  page  329]: 


(3.28) 


where  Fdrag=  drag  force  (N),  p—  atmospheric  density  (kg/rn3),  V=  velocity  (m/s), 
Cp=  drag  coefficient  and  Ap=  projected  area  (nr2). 

Since  the  satellite  is  assumed  to  be  in  a  circular  orbit,  the  tangential  velocity 
V  is  constant  which  is  given  by  [11,  page  70]: 


(3.29) 


where  the  standard  gravitational  parameter  p=  398600^-  and  the  orbital  radius 


v orbit =  6378.135/cm  +  Altitude.  Note,  the  velocity  vector  is  in  the  — Si  direction. 


The  density  p  depends  on  altitude  and  is  calculated  by  the  Matlab®  code  in 


Appendix  C.  From  Moe’s  research,  a  good  experimental  approximation  for  Cp  is 
2.2  [4,  page  4],  The  projected  areas  Ap  and  their  locations  with  respect  to  the  CG 
are  the  only  parameters  that  can  be  controlled  after  the  spacecraft  is  in  orbit. 

Altitudes  above  125  km  are  in  the  free  molecular  flow  regime  [1,  page  316]. 
In  the  free  molecular  flow  regime,  particles  are  typically  modeled  either  as  specular 
or  diffuse  reflections.  A  specular  reflection  assumes  that  molecules  are  perfectly 
clastic  where  the  tangential  velocity  is  constant  and  the  normal  velocity  is  equal  and 
opposite  before  and  after  reflection.  The  diffuse  model  assumes  the  molecules  are 
reflected  in  a  diffuse  manner  and  have  no  memory  of  previous  velocities.  Either  model 
imparts  a  force  normal  to  the  surface  and  the  reflections  are  wrapped  up  into  the 
drag  coefficient  which  as  stated  earlier  was  determined  experimentally.  Spacecraft 
in  LEO  always  experience  a  small  drag  torque  since  it  is  not  possible  to  construct  a 
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spacecraft  that  has  its  CG  at  the  exact  geometric  center.  Some  assumptions  made 
for  the  following  drag  torque  equations  are: 

•  The  drag  panels  are  visible  to  the  incoming  molecules  on  the  front/outward 
facing  surfaces. 

•  The  drag  coefficient  was  assumed  to  be  constant  at  2.2  [4,  page  4], 

•  Atmospheric  density  is  constant  at  a  constant  altitude  and  averaged  over  the 
orbit. (neglecting  solar  effects  and  atmospheric  perturbations). 

•  The  incoming  molecular  velocity  is  equal  and  opposite  the  orbital  velocity 
(non- rotating  atmosphere) . 

•  The  mass  of  the  control  panels  and  arms  are  negligible  compared  to  the  mass 
of  the  satellite. 

3-4-1  Spacecraft  Configuration.  The  linearized  model  consists  of  a  basic 
cube  shaped  satellite  bus  with  drag  panels  that  can  be  extended  into  the  airstream 
for  pitch  and  yaw  control.  Two  wedges  acting  as  “ailerons” ,  to  produce  torques  about 
the  roll  axis  (See  Figure  2.9).  Wedges  were  chosen  instead  of  flat  panels  to  simplify 
the  linearization  of  the  equations  of  motion.  When  linearizing  the  torque  equations 
about  zero  degrees  for  the  ailerons,  all  variables  would  become  higher  order  terms 
since  the  projected  areas  and  linearization  angles  were  zero.  Since  Higher  order  terms 
were  to  be  neglected,  there  were  no  variables  left  in  the  equations.  By  modelling  the 
spacecraft  with  wedge  shaped  ailerons,  the  linearization  problem  goes  away  as  does 
the  heaviside  function  in  projected  area  equations. 

3-4-2  Drag  Effects  from  Spacecraft  Body.  The  spacecraft  body  doesn’t  pro¬ 
duce  any  torques  due  to  drag  unless  the  CG  is  not  at  the  geometric  center  (assuming 
a  symmetric  spacecraft).  Although,  in  the  linearized  torque  equations,  the  CG  loca¬ 
tion  is  taken  to  be  at  the  geometric  center  thus  drag  effects  from  the  spacecraft  body 
are  neglected  in  the  linearized  model.  In  the  nonlinear  equations  in  Appendix  F,  the 
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CG  location  can  be  moved  off-center  to  produce  small  torques.  No  further  analysis 
was  done  with  drag  from  the  spacecraft  body.  The  remaining  analysis  focused  on 
the  control  panels. 


3-4-3  Drag  Effects  from  Spacecraft  Drag  Panels.  The  first  step  was  to  derive 
the  rotation  matrix  to  transform  the  drag  panel  reference  frames  (the  F)  frame)  to 
the  body  frame.  Each  drag  panel  rotation  was  based  on  equation  3.30 
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1 _ 

s(h) 
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0 

to 

0 

0 

1 

(3.30) 


Each  drag  panel  has  only  one  degree  of  freedom,  therefore  equation  3.30  is  simplified 
requiring  one  rotation  per  panel.  Each  panel’s  rotation  matrix  contains  a  single 
rotation  about  its  respective  axis.  Figure  3.3  shows  how  each  arm  rotates  with 
respect  to  the  satellite  bus,  i  is  either  the  top,  bottom,  left,  or  right  arm.  Equation 
3.30  for  each  panel  simplifies  to 


Rfb 


Top 


c(  (j)  Top  Arm)  0  s(  (bj1  op  Arm) 
0  1  0 

)  0  c(  f  Top  Arm) 


(3.31a) 


r>FB  _ 
a  Bot  — 
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Cyf'BotArm)  0  ^i^fBotArmf) 

0  1  0 

S^fBotArm)  0  c(^(f>BotArm) 

c(0  Left  Arm)  s((j)LeftArrn)  0 

Left  Arm)  c(0  Left  Arm)  0 

0  0  1 


(3.31b) 


(3.31c) 
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rfb 


Right 


c(  (pRightArm)  s(  (pRightArm )  0 

s(  ( pRightArm )  c(  4*  Right  Arm)  0 

0  0  1 


(3.31(1) 


For  the  wedges,  the  angles  of  the  faces  are  found  by  adding  or  subtracting  the  wedge 


Figure  3.3:  Angles  for  each  control  arm  with  respect  to  body 
frame 


angle  to  get  the  total  angle.  The  angles  are  defined  as  follows: 


<5 h  eftAUTop  (pLeftAil  ^pWedge  (3.32a) 

(pLeftAUBot.  =  (pLeftAil  +  (pWedge  (3.32b) 

( pRightAUTop  (pRightAil  (pWedge  (3.32c) 

(p Right AilBot  0  RtghlAil  T  (pWedge  (3.32d) 


where  the  aileron  control  angles  are  (pmghtAu—  ~(pLeftAU ,  and  (pWedge  is  half  the  wedge 
angle.  See  Figure  3.4.  The  rotations  for  each  surface  then  become 
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tdFB  _ 
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(3.33d) 

where  the  control  arm  angles  O  Arm  are  Ot  op  Arm  •  0  Boi.Arm-  (pLeftArmi  and  (j)  Right  Arm  i 
are  the  top,  bottom,  left  and  right  arms  that  are  connected  to  the  drag  panels.  The 
angles  between  the  b\  direction  and  the  arms  range  from  0  —  90°.  In  equations  3.31a 
and  3.31d,  the  arm  angles  have  a  negative  sign  so  that  all  angles  will  be  positive  to 
make  it  easier  to  remember  that  all  angles  to  be  input  to  the  controller  are  positive. 
Two  of  the  rotation  angles  in  each  equation  are  always  zero  since  there  is  only  one 
degree  of  freedom  for  each  panel. 


The  rotations  from  the  A  frame  to  the  F  frame  are  used  to  determine  normal 
force  for  each  panel  due  to  the  drag  in  the  A  frame.  This  is  found  by  RtA  =  RFBRBA 
and  thus  becomes: 
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(3.34g) 
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RFA  RightAilBot  ~  RFB Right AilBotRBA  (3.34h) 

The  angles  between  the  incoming  molecules  and  the  inward  normal  directions 

are: 
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The  projected  areas  from  the  panels  become: 


ATop  =  L2 H (cos(aTop))  cos (aTop)  (3.36a) 

ABot  =  L2 H (cos(aBot))  cos (aBot)  (3.36b) 

ALeft  =  L2 H (cos(aLeft))  cos(aLeft )  (3.36c) 

Ajught  R  H (cos(cr Right)')  cos(o:/jjg/jj)  (3.36d) 
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where  L  is  the  side  length  of  each  square  drag  panel,  Lau  is  the  length  (span)  of 
each  aileron  and  Wau  is  the  width  of  each  aileron.  Calculate  forces  from  the  drag 
panels  in  the  A  frame: 


Ftopa  =  2 CdAToVpV 2v 
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where  v  =  is  the  velocity  unit  vector  and  is  defined  as: 
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Calculate  force  from  drag  panels  in  the  F  frame  and  pull  out  the  inward  normal 
component: 
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Convert  the  inward  normal  component  of  each  surface  in  the  F  frame  to  the  B  frame: 


FTopb  —  {RFBToP)T  Fpopp  (3.40a) 

FBotB  =  {RFB Bot)T  FsotF  (3.40b) 
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Find  r  from  the  CG  to  the  center  of  pressure  for  each  panel.  (See  Figure  3.4 
for  aileron  details.)  (Note:  CG  is  at  the  geometric  center  for  the  linearized  model, 
i.e.  CG=0): 


V  Top 


TBot  = 


rLeft  ~ 


^ Right 


0-5 L  (L Arm  0.5 L panel)  COS  Arm) 

0 

0.5L  (i,4rm  F  0. 5 L panei)  sill ((f^TopArm) 

0.5L  [F  Arm  +  0.5 L Panel)  COS {(f>BotArm) 

0 

0-5L  +  [F  Arm  ~h  0. 5 L panel)  sill(0 BotArm ) 

0.5L  [F  Arm  +  0.5 Lpanel)  COS [(j)  Left  Arm) 
0.5L  (i'Arm  “1“  0.5Lpane/)  sill(0  Left  Arm) 

0 

0.5L  (-^Am  H-  0. 5 L Panel)  COS (0 Right Arm) 

0.5  L  +  {Fa  rm  4“  0.5  L panel)  si  Right  Arm) 

0 


CG 


- CG 


- CG 


CG 


(3.41a) 


(3.41b) 


(3.41c) 


(3.41d) 
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^  LeftAilTop 


-CG 


(3.41e) 


0.5  Wau  sill  (pWedge  sill  4*  Left  Ail 
—  0.5  (Lside  +  Lau) 

0 .5W Ail  sill  (f>Wedge  COS  (pLeftAil 


L  LeftAilBot 


0.511  4 ,1  sill  (pWedge  sin  (pLeftAil 
0.5(Z/1g,jc;g  -|-  Lau) 

0-51 1  ~\il  sill  4>Wedge  COS  (j) Left  Ail 


-CG 


(3.41f) 


^ Right  AilTop 


O.BWau  sill  (pWedge  sill  4  Right  Ail 
0.5(Z/gj^e  -|-  Lau ) 

0 . 5  M,4?7  sill  (pWedge  COS  0  flight  Ail 


-CG 


(3.41g) 


T RightAilBot 


0.5W^4j/  sill  0  Wedge  sill  $  Right  Ail 


0-5  (LSide  +  Lau) 


-CG 


(3.41h) 


0.511  \  \il  sill  0W  ed.ge  COS  (j)  Right  Ail 


where  Lside  is  the  side  length  of  the  spacecraft  body,  Lat m  is  the  length  of  the 


LOCATION  OF  THE 
CENTER  OF  PRESSURE 
ALONG  THE  bl  AND  b3 
DIRECTIONS 
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control  arm,  Lpanei  is  the  side  length  of  the  square  panels,  (pLArm ,  (pRArm ,  (pBArm  and 
(pTArm  are  the  angles  of  the  arms  from  their  retracted  positions,  (p  Left  Ail  and  (p  Right  Ail 
are  the  angles  of  the  ailerons,  and  (pwedge  is  half  the  angle  of  the  wedge. 

Finally,  the  torque  equations  for  the  drag  panels  with  all  the  substitutions 
become: 


FI Top  ^ T  opXF TopQ 

(3.42a) 

M-Bot  ^ BotX-FBotB 

(3.42b) 

FI  Left  T*  Left^-F^eft  b 

(3.42c) 

FI  Right  ^  Right^-F Right  b 

(3.42d) 

FI  Left  AilT op  ^ Left  AilT op^-F Left  AilT  op b 

(3.42e) 

FI  LeftAilBot  ^  Left  AilBot%~F Left  AilBot  b 

(3.42f) 

FI Right  AilT op  ^ Ri g  ht  AilT  op^-F Right  AilT  op  b 

(3.42g) 

FlRightAilBot  T Right  AilBot^F Right AilBot  r 

(3.42h) 

To  get  the  total  torques  for  the  drag  control  panels,  equations  3.42  a-h)  are  added: 

M Panels  =  Mpop  +  Mpot  +  Mpeft  +  MRight  +  MpeftAilTop 

T  M LeftAilBot  T  Ld Right AilT op  4“  dl Right AilBot  (3.43) 

If  the  torques  from  the  body  were  included,  they  would  be  added  here,  however  in 
this  case,  MSatBody  =  0: 


Mpotal  —  MsatBody  +  Mpaneis  (3.44) 

3-4-4  Linearized  Torques  from  Control  Panels.  Linear  analysis  required 
the  nonlinear  torque  equations  to  be  linearized.  By  having  the  ailerons  as  wedges 


33 


and  assuming  small  angles,  the  outward  facing  control  surfaces  were  always  facing 
the  incoming  molecules.  Thus,  the  heaviside  function  was  no  longer  necessary  for 
the  linear  model.  Each  control  arm  for  the  pitch  and  yaw  controls  were  linearized 
about  (p oArrn,  which  is  the  linearization  angle  of  the  control  arms.  The  ailerons  were 
linearized  about  <poAU,  where  4> oAil  =  0°  (aligned  with  the  b\  axis).  After  using 
both  MathCad  and  Mat  lab®  symbolic  solvers,  the  combined  and  linearized  torque 
equations  were  found: 


Mi 

CRAilS(j)RightAil  +  ClAU^  Left Ail) 

m2 

= 

C BArm^&BotArm  T  CR  Arm^  0T  op  Arm  T  Cq2@2 

m3_ 

C LArm^^  Left  Arm  T  T  Cq39^ 

where  CRAu,  CLAu,  CTArm ,  CBArm,  CLArm,  CRArm ,  Ce2,  and  Cq 3  are  dehned  as  follows: 


CRAu  =  [0.5 LAu2WAu  Sm((/)wedge) 

T  ().5LAii W AH Lh2  sin ((ftw  edge) 

1 .5LAh~WAh  edge)  COS  (</> Wedge) 

-  l.5LAaWAULb2sm((j)Wedge)cos2((l)Wedge)]CDpV2  (3.46a) 

Clau  =  —CRAu  (3.46b) 


CTArm 


[0.5 Lpanei  W panel  COS{(f)0Arm)  sin(0OArm) 

C panel  ^^panel  L Arm  COs(0oArm)  sin(0oArm) 

0.75  Lpanel  W panel  Lbl  COS  (00  Arm)  sin(0o  Arm) 
0.75 Lpanei W panei Z/&3  COS  (00 Arm) 

0  •  75Lpanei  W panel  Lb3  COs(0o^rm) 

0.25 Lpanel^^ panel l^bi  sin(0o Arm)\^DpV 


(3.46c) 
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—  Cl  Arm 


(3.46c!) 


Cl  Arm  [  0.5Z/pane/  ^ panel  COs(0q  Arm)  sin(0OArm) 

L panelW panel L Arm  COs(0o Arm)  S1h(0o Arm) 

0.75 LparielW panel Lbx  COS  (00 Arm)  Slfr(0OArm) 
~h  0.75Z/pane/ WpanelLb2  COS  (00 Arm) 

0 .75Lpanel\VpanelLl)2  COS (00 Arm) 

~b  0 •^7)Lpanei\VpaneiL})1  sin(0oAr,m )] CdpV 


CftArm  —  C LArm 


Cq2  [2 LaH^V Ail  COS  Wedge)  sill(0ip^ edge) 

2LAuWau  COS ((j)W edge)  sill(0p^e^e) 

Lpanel  2Wpanel  COs(0oArm)  Sm(0oArm) 

2 Lpanei  W panei  LArm  COs(0oArm)  sin(0oArm) 
CpanelW panel Cbi  COS  (00 Arm)  sill(0oArm) 

“1“  Cpane{\VpanelL}D 3  COS  (00 Arm) 

LpanelW^ panel  Lbs  COs(0o  Arm)\CDpC 


Cq-s  [Lpanel  WpanelLb2  COS  (0OArm) 

2  LpaneiWpanei  LArm  COs((/>0 Arm)  S^1(<i^0Arm) 
LpaneiWpanei Lbi  COS^(0oArm)  sin(0oA7,m) 

Lpanel  2Wpanel  COs(0OApm)  sin (00 Arm) 

LpanelW panel Lb$  COs(0o^rm )] Cj^pV 


(3.46e) 

(3.46f) 


(3.46g) 


(3.46h) 
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Lau  is  the  spanwise  length  of  the  aileron  surface,  Wau  is  the  chord  length  of  the 
aileron  surfaces,  Lpanei  is  the  length  of  the  drag  panels,  Wpanei  is  the  width  of  the 
drag  panels,  Lbl  is  the  length  of  the  satellite  body  in  the  bi  direction,  Lb2,  is  the 
length  of  the  satellite  body  in  the  b 2  direction,  and  Lbs,  is  the  length  of  the  satellite 
body  in  the  63  direction. 


3.5  Combined  Linearized  Equations  of  motion 

The  final  step  in  deriving  the  EOM  was  to  combine  all  the  equations.  By 
plugging  the  linearized  torque  equations  3.45  into  the  overall  equations  of  motion 
3.27  and  simplifying,  the  equations  become: 


01 

01 

01 

h 

-  C'-Y,  ■ 

02 

+  Cx2  ■ 

02 

A 

A 

A 

b0T  opAr  m 
80BotArm 
8  0  Left  Arm 
80RightArm 

80LeftAil 

80RightAil 


(3.47) 


where  80i  is  the  deflection  angles  for  the  respective  control  arms  and  ailerons; 
80TopArmi  80BotArrm  80LeftArmi  nighl.Arrn-  8  0  [.eft  Ail  •  8  0  Right  Ail  i  and  C X\  ,  C \2 ,  and 
Cx3  are  defined  as  follows: 
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(i C  -  B)n2 0 

0 

0 

1 

B 

0 

0 

P 

to 

O 

0 

0 

1 

c_ 

0 
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(3.48) 


(3.49) 
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1 

A 

0 

0 

0 

0 

0 

0 

Cl  Ail 

1 — 

Cj 

Cx3  - 

0 

1 

B 

0 

CxArm 

C BArm 

0 

0 

0 

0 

0 

0 

1 

c 

0 

0 

C^LArm 

C R Arm 

0 

0 

Equation  3.47  with  all  the  substitutions  is  the  final  equation  to  be  modeled  in 
Simulink®. 


37 


IV.  Dynamic  Model 


4.1  Dynamic  Model  Overview 

The  linearized  dynamic  equations  3.47  were  modeled  in  Simulink®.  All  terms 
were  left  as  variables.  Inputs  so  different  scenarios  could  be  easily  simulated  by 
changing  the  input  variables.  If  more  details  of  a  satellite  are  known,  such  as  it’s 
dimensions,  mass  and  orbital  altitude,  this  model  can  be  easily  updated  to  simulate 
a  variety  of  satellites.  The  only  restriction  is  that  the  satellite  has  to  have  the  same 
basic  geometric  design. 

4-2  Dynamic  Model  Validation 

The  following  is  the  validation  of  the  open  loop  and  closed  loop  model  by  using 
eigenvalues  and  eigenvectors,  where  the  eigenvalues  are  the  frequencies  (cycles/orbit) 
of  this  system  and  the  eigenvectors  determine  initial  conditions  that  will  predict 
system  response.  The  results  of  this  analysis  show  how  the  general  shape  of  the 
response.  Such  responses  could  be  pure  oscillatory  or  exponential.  In  retrospect, 
the  eigenvectors  could  have  been  scaled  so  the  amplitude  of  the  responses  remained 
small,  however,  having  a  larger  amplitude  response  does  not  alter  the  outcome  of 
this  section.  All  simulations  performed  in  this  section  were  run  at  300km  altitude 
and  the  density  was  set  to  zero  for  the  no  torque  case  and  the  density  was  calculated 
for  300km  altitude  for  the  case  with  torques. 

4-2.1  Dynamic  Model  Validation  Without  Torques.  The  first  step  was  to 
set  all  external  torques  of  equation  3.27  to  zero.  The  rest  of  equation  3.27  was 
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rewritten  in  a  state  space  ( X  =  A  ■  X)  format  as  follows. 


01 
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0i 

e2 
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<^b. 
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0 

0 

0 

0 

0 

1 

03 
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C+A-Bn 
A  n 

0i 
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0 

0 

0 

0 

0 

02 

o. 

0 

0 

A~Bn 2 

B~A~cn 

0 

0 

1 

CO 

*^b 

(4.1) 


From  the  state  space  equation  (equation  4.2),  the  pitch  axis  (02)  has  no  influence  on 
the  other  axes,  hence  it  is  decoupled  from  the  roll  (6*i)  and  yaw  (03)  axes.  02  was 
shown  to  be  stable  based  on  test  simulations.  Since  02  is  stable  and  decoupled,  all 
references  to  d2  were  removed  leaving  the  following  equation. 


0i 
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0i 

CO 

0 
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^b 

CO 

0i 

C—B  2 

A  H 
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0 

C+A-B 

A  n 
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.  co 
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0 

A~Bn 2 

B~A~cn 
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1 

CO 

•Qb 

(4.2) 


The  next  step  was  to  test  the  model  with  different  combinations  of  major,  minor 
and  intermediate  axes  for  the  principal  moments  of  inertia.  If  the  spacecraft  is  nadir 
pointing,  the  pitch  axis  will  always  be  rotating  with  the  orbit.  If  the  pitch  axis  is 
either  the  minor  or  major  axis  of  rotation,  the  spacecraft  will  be  stable.  If  the  pitch 
axis  is  the  minor  axis,  the  spacecraft  will  be  in  an  unstable  configuration.  Several 
combinations  of  MOI  were  simulated  where  the  moments  of  inertia  about  the  roll, 
pitch  and  yaw  axes  are  A,  B,  and  C  respectively.  For  the  first  simulation,  the  b2 
axis  was  set  as  the  minor  axis,  where  A  =  20,  B  =  10,  and  C  =  30.  The  calculated 
eigenvalues  A,  are  arranged  in  a  diagonal  matrix  of  eigenvalues  (A)  and  the  calculated 
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eigenvectors  are  the  columns  of  (e). 


A  = 


i  ■  n 

0 

0 

0 


0  0 

-i  ■  n  0 

0  i  ■  n 
0  0 


0 

0 

0 

-i  ■  n 


(4.3) 


e  = 


—i 

i 

-0.865* 

0.865* 

1 

1 

1 

1 

n 

n 

0.500n 

0.500n 

i  ■  n  -i-n 

0.576n 

— 0.576r* 

(4,4) 


The  initial  conditions  ( 9i  and  9i )  for  the  first  and  second  complex  conjugate 
pair  of  eigenvectors  can  be  found  by: 
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1.730 

2 

n 

-1.156n 


(4.5) 


(4.6) 


When  using  the  initial  conditions  from  equation  4.5  in  the  Matlab®  code  in  Appendix 
B,  the  response  for  each  eigenvalue  and  associated  eigenvectors  should  be  purely 
oscillatory  (since  it  has  only  imaginary  components)  and  have  a  period  of  1  cycle 
per  orbit. 
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Figures  4.1,  4.2  and  4.3,  shows  the  angular  position,  angular  rate  and  angular 
acceleration  of  this  system.  The  system  response  shows  that  the  oscillations  are 
purely  sinusoidal  with  a  period  equal  to  the  orbital  period  of  5431  seconds.  Recall 


that  the  roll,  pitch  and  yaw  axes  are  6h,  62,  and  63  respectively. 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  4.1:  Angular  Position  or  first  pair  of  eigenvectors  with 
A  =  20, B  =  10,  and  C  =  30 
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Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  4.2:  Angular  Rate  for  first  pair  of  eigenvectors  with 

A  =  20,  B  =  10,  and  C  =  30 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  4.3:  Angular  Acceleration  for  first  pair  of  eigenvectors 
with  A  =  20,  B  =  10,  and  C  =  30 
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For  the  second  set  of  complex  conjugate  eigenvectors,  there  should  be  0.576 
cycles  per  orbital  period.  The  period  of  oscillations  is  therefore  9429  seconds  which 
can  be  seen  from  Figures  4.4,  4.5  and  4.6. 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  4.4:  Angular  Position  for  second  pair  of  eigenvectors 

with  A  =  20,  B  =  10,  and  C  =  30 
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2 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


x  Id"3 


Figure  4.5:  Angular  Rate  for  second  pair  of  eigenvectors  with 
A  =  20,  B  =  10,  and  C  =  30 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  4.6:  Angular  Acceleration  for  second  pair  of  eigenvec¬ 
tors  with  A  =  20,  B  =  10,  and  C  =  30 
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Another  run  was  done  with  62  as  the  intermediate  axis,  where  A  =  10,  B  —  20, 
and  C  =  30.  In  this  configuration,  the  satellite  will  be  unstable.  This  was  tested 
to  evaluate  different  combinations  of  the  MOI  to  test  different  possibilities.  The 
eigenvalues  and  eigenvectors  for  this  particular  MOI  are: 
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(4.7) 


(4.8) 


Once  again,  the  same  procedure  was  done  to  obtain  the  initial  conditions  for  this  set 
of  eigenvalues  and  eigenvectors.  The  initial  conditions  become 
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(4.10) 


The  case  where  the  initial  conditions  from  equation  4.9  are  used,  the  outputs  for 
6\  and  d3  in  position,  angular  rate  and  angular  acceleration  should  be  sinusoidal 
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oscillations  at  one  cycle  per  orbit  due  to  the  eigenvectors  containing  only  imaginary 
components.  The  sinusoidal  oscillations  can  be  seen  in  Figures  4.7,  4.8,  and  4.9 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  4.7:  Angular  Position  for  first  pair  of  eigenvectors  with 
A=10,B  =  20,  and  C  =  30 
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Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  4.8:  Angular  Rate  for  first  pair  of  eigenvectors  with 

A  =  10,  B  =  20,  and  C  =  30 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  4.9:  Angular  Acceleration  for  first  pair  of  eigenvectors 
with  A  —  10,  B  —  20,  and  C  =  30 
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The  second  set  of  eigenvalues  and  eigenvectors  have  no  imaginary  parts,  there¬ 
fore  the  0’s  will  grow  exponentially.  Exponential  growth  for  the  angles  is  due  to  the 
spacecraft  being  in  an  unstable  configuration  since  the  pitch  axis  is  the  intermediate 
MOI  axis.  The  results  of  this  run  can  be  seen  in  Figures  4.10,  4.11,  and  4.12. 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  4.10:  Angular  Position  for  second  pair  of  eigenvectors 
with  A  —  10,  B  —  20,  and  C  =  30 
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Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  4.11:  Angular  Rate  for  second  pair  of  eigenvectors  with 
A  =  10,  B  =  20,  and  C  =  30 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  4.12:  Angular  Acceleration  for  second  pair  of  eigenvec¬ 
tors  with  A  —  10,  B  —  20,  and  C  =  30 
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4-2.2  Dynamic  Model  Validation  With  Torques.  The  next  step  was  to  add 
the  atmospheric  density  to  the  model  and  perform  a  similar  eigenvalue/eigenvector 
analysis  to  see  how  the  model  responds.  The  same  procedure  was  carried  out  as  in 
the  previous  section.  Equation  3.47  including  in  state  space  format  ( X  =  A-x  +  B-u ) 
including  torques  is  as  follows: 
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(4.11) 


Once  again,  since  62  is  totally  decoupled  from  9 1,  03,  9\  and  6s  all  references  to 
62  were  removed  leaving  the  following  A  matrix  to  do  the  eigenvalue/eigenvector 
analysis. 
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The  eigenvalues  and  eigenvectors  from  Equation  4.12  were  calculated  using  the 
Matlab®  code  in  Appendix  B.  This  simulation  was  also  done  at  300km  altitude 
however,  atmospheric  density  was  p  =  2.4  x  10~n  Note,  the  most  stable  satel¬ 
lite  configuration  is  with  62  as  the  major  MOI  axis,  63  as  the  minor  axis  and  b±  as 
the  intermediate  axis.  The  following  analysis  was  done  with  A  =  20,  B  =  30,  and 
C  =  10. 
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(4.14) 


Once  again,  by  taking  each  complex  conjugate  pair  of  eigenvectors,  the  initial  con¬ 
ditions  can  be  found. 
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Upon  running  the  Simulink®  model  with  the  initial  conditions  computed  from  equa¬ 
tion  4.15,  the  results  indicate  that  the  model  is  correct  since  the  frequency  of  oscil¬ 
lation  should  again  be  one  cycle  per  orbit.  See  Figures  4.13,  4.14,  and  4.15 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  4.13:  Angular  Position  with  torque  for  first  pair  of 

eigenvectors  with  A  =  20,  B  =  30,  and  C  —  10 
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Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  4.14:  Angular  Rate  with  torque  for  first  pair  of  eigen¬ 
vectors  with  A  =  20,  B  =  30,  and  C  —  10 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  4.15:  Angular  Acceleration  with  torque  for  first  pair 

of  eigenvectors  with  A  =  20,  B  =  30,  and  C  —  10 
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The  same  analysis  was  done  on  the  remaining  set  of  initial  conditions  from 
equation  4.16.  Figures  4.16,  4.17,  and  4.18,  show  that  the  cyclic  period  of  oscillations 
is  reduced  (higher  frequency).  Since  the  period  is  approximately  315s,  the  frequency 
can  be  calculated  from  uj  =  djT  =  0.0199^,  which  is  close  enough  to  the  eigenvalue 
0.0253  for  model  validation. 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  4.16:  Angular  Position  with  torque  for  second  pair  of 
eigenvectors  with  A  =  20,  B  =  30,  and  C  =  10 
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Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  4.17:  Angular  Rate  with  torque  for  second  pair  of 

eigenvectors  with  A  =  20,  B  =  30,  and  C  —  10 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  4.18:  Angular  Acceleration  with  torque  for  second  pair 
of  eigenvectors  with  A  =  20,  B  =  30,  and  C  —  10 
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4-3  Proportional  Derivative  Controller 

Since  the  Simulink®  model  is  linearized,  the  control  deflections  are  kept  rea¬ 
sonably  small.  The  controller  is  a  proportional  controller  with  the  constants  selected 
in  such  a  way  that  the  control  deflections  remain  small  due  to  the  linearized  torque 
equations.  As  discussed  in  chapter  V,  the  control  deflections  do  in  fact  remain  small 
thus  giving  a  more  realistic  result.  With  the  controller,  the  roll,  pitch  and  yaw  axes 
will  eventually  settle  down  in  the  correct  satellite  orientation.  Originally,  the  con¬ 
troller  was  set  up  such  that  the  control  arm  outputs  were  equal  to  some  constant 
multiplying  the  angular  rates.  The  roll  axis  seemed  to  take  a  while  to  return  to 
zero  degrees.  The  only  reason  the  roll  axis  would  return  to  zero  was  because  it  was 
coupled  with  the  yaw  axis  and  the  dynamics  of  the  system  forced  both  the  roll  and 
yaw  axis  to  work  together  to  align  the  body  frame  to  the  orbit  frame.  Because  of 
the  long  settling  time,  the  controller  was  modified  so  that  each  control  arm  inputs 
for  each  axis  was  equal  to  a  constant  multiplying  the  angular  rates  plus  another 
constant  multiplying  the  angle  between  the  body  and  orbit  reference  frames  (see 
equation  4.17.  Adding  the  constant  times  the  anglular  position  brought  the  roll  axis 
back  to  zero  degrees  in  a  more  reasonable  time  frame. 

5(f)  =  C9  +  0.19  (4.17) 
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V.  Results 


5.1  Altitude  Variations 

Based  on  equation  3.28,  drag  is  proportional  to  the  square  of  the  velocity  and 
is  proportional  to  the  atmospheric  density.  A  spacecraft  in  LEO  will  experience 
orbital  decay  due  to  drag  which  will  lower  its  orbit  over  time.  As  orbital  altitude  is 
reduced,  both  the  density  and  velocity  increase  thus  having  a  significant  increase  on 
the  drag  force.  Therefore,  the  most  significant  variable  that  affects  drag  is  altitude. 
For  the  following  subsections,  the  simulations  were  run  at  altitudes  varying  from 
200km  to  600km  in  100km  increments.  At  these  altitudes,  the  density,  velocity  and 
orbital  rate  were  calculated  and  are  shown  in  Table  5.1.  Table  5.2  shows  all  the 
inputs  to  the  model  for  all  the  simulations  where  altitude  is  varied.  The  MOI  used 
for  this  section  are  shown  in  equation  5.1  (See  Appendix  D.3  for  details  on  how  this 
was  determined).  Table  5.2  shows  the  inputs  to  the  Simulink®  model.  Note,  only  a 
few  plots  were  included  in  this  chapter,  however,  all  the  plots  have  been  included  in 
Appendix  E. 


Table  5.1:  Density,  velocity  and  orbital  rate  calculated  from  Simulink®. 


Altitude  (km) 

Calculated  Density  (®tt) 

Velocity  (2) 

Orbital  Rate  (^) 

200 

2.79xl0~lu 

7784 

0.001183 

300 

2.42xl0~n 

7726 

0.001157 

400 

3.73xl0”12 

7669 

0.001131 

500 

6.97xl0~13 

7613 

0.001107 

600 

1.45xl0~13 

7558 

0.001083 

MOI  = 


A 

0 

0 

92 

0 

0 

0 

B 

0 

= 

0 

107 

0 

0 

0 

c 

0 

0 

68 

(5.1) 


The  proportional  derivative  controller  constants  multiplying  the  angular  rates  are  50 
for  roll,  20  for  pitch  and  20  for  yaw.  The  constants  multiplying  the  angle  between  the 
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Table  5.2:  Variables  input  into  the  Simulink®  model  for  all  the  altitude  variation 
simulations. 


Variable 

Value 

L  Panel 

lm 

^Vpanel 

Irri 

L  Arm 

1  m 

45o 

Lau 

lm 

Wau 

lm 

lm 

-^62 

lm 

T&3 

lm 

4*  W  edge 

22.5° 

Mass 

500% 

Cd 

2.2 

body  and  orbit  frames  are  0.1.  For  the  initial  conditions,  the  angular  displacement 
of  the  satellite  body  with  respect  to  the  orbital  frame  are  chosen  to  be  offset  10° 
and  the  angular  rates  set  to  The  linearization  angle  for  the  control  arms  was 

chosen  to  be  45°.  Several  factors  were  considered  in  selecting  the  angle.  First,  when 
the  panels  are  extended,  only  the  windward  facing  sides  are  subject  to  being  hit 
by  molecules.  Since  small  angles  are  assumed,  the  back  side  of  the  panels  will  not 
be  exposed  to  the  incoming  molecules.  The  small  angle  assumption  between  the 
orbit  and  body  frames  simplified  the  linearization  since  the  heaviside  function  was 
no  longer  needed.  Another  reason  is  that  at  smaller  linearization  angles,  the  wedge 
shaped  ailerons  would  shadow  the  left  and  right  panels  used  for  yaw  control.  The 
ailerons  for  the  nonlinear  case  would  not  have  this  problem  since  the  ailerons  would 
consist  of  flat  panels,  not  wedges. 

5.1.1  Simulation  at  200  km  Altitude.  For  this  simulation,  the  inputs  to 
the  modeled  are  represented  in  Table  5.2  with  an  altitude  of  200 km.  The  density, 
orbital  velocity  and  orbital  rate  for  this  altitude  are  summarized  in  Table  5.1.  The 
first  simulation  only  had  the  roll  axis  displaced,  the  2nd  was  pitch  and  the  3rd  was 
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the  yaw  axis.  Then  a  fourth  simulation  was  run  where  all  three  axes  of  roll  pitch 
and  yaw  were  displaced  simultaneously.  By  comparing  Figures  5.1,  and  5.4,  it  can 
be  seen  that  roll  and  yaw  are  coupled  and  do  have  some  affect  on  each  other.  When 
both  yaw  and  roll  are  offset  simultaneously  as  shown  in  Figure  5.6,  the  coupling  is 
more  pronounced. 

5. 1.1.1  200km  Altitude,  Roll  Offset  by  10°.  The  first  simulation  was 

done  where  roll  was  offset  by  10°.  To  compare  the  settling  times  for  the  various 
altitudes,  the  time  constant  must  be  found.  The  time  constant  r  is  the  time  it  takes 
for  the  amplitude  to  decay  by  A"ff° ,  where,  Amp o  is  the  initial  amplitude  or  initial  9, 
and  e  ~  2.71.  The  time  constant  can  be  approximated  directly  from  the  plots  once 
the  initial  amplitude  is  known.  For  the  simulations,  the  initial  conditions  determine 
Ampo,  and  thus  Ampo  =  10°  or  0.175rad.  The  amplitude  for  one  time  constant  is 

Anrpo  =  Q.064. 
e 

Roll,  Pitch,  &  Yaw  Position  vs.  Time  Control  Angle  Deflections 


(a)  (b) 

Figure  5.1:  (a)  Angular  displacement  with  roll  offset  10°. 

(b)  Control  angles  with  roll  10°. 

Since  the  roll  doesn’t  oscillate,  the  time  constant  can  be  read  from  the  plot 
where  r  500s 
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5. 1.1.2  200km  Altitude,  Pitch  Offset  by  10°.  The  pitch  is  decoupled 
from  the  roll  and  yaw  axes  therefore  the  only  oscillations  should  be  in  the  pitch  axis. 
Figure  5.2(a)  shows  this  to  be  true. 


(a)  (b) 

Figure  5.2:  (a)  Angular  displacement  with  pitch  offset  10°. 

(b)  Control  angles  with  pitch  10°. 


The  next  step  was  to  pick  points  on  the  plot,  and  fit  them  to  a  curve.  The  time 
constant  for  the  best  fit  can  then  be  determined.  Since  the  pitch  axis  oscillates,  the 
locations  of  the  peaks  were  used  for  the  curve  fit  algorithm  in  Appendix  B.  Assuming 
an  exponential  decay,  the  plot  should  fit  the  following  equation: 

Y  =  [exp(—^-)\  *a  (5.2) 

r 

The  plot  shown  in  Figure  5.3  closely  approximates  the  data  points  and  the  time 
constant  is  shown  to  be  ~  130s 
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Figure  5.3:  Curve  fit  for  200km  altitude  with  pitch  offset  10° 


5. 1.1.3  200km  Altitude,  Yaw  Offset  by  10°.  The  same  procedure 
performed  in  the  previous  section  was  applied  to  the  simulation  where  the  yaw  was 
offset  by  10°.  The  results  are  shown  in  Figure  5.4  (a)  and  5.5.  Here,  the  time 
constant  is  ~  85s. 
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Control  Angle  Deflections 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


(a) 
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Figure  5.4:  (a)  Angular  displacement  with  yaw  offset  10°. 

(b)  Control  angles  with  yaw  10°. 


Figure  5.5:  Curve  fit  for  200km  altitude  with  yaw  offset  10° 
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5. 1.1. 4  200km  Altitude,  Roll,  Pitch  and  Yaw  Offset  by  10°.  For  this 

simulation,  all  3  axes  were  offset  simultaneously  which  would  be  a  more  realistic 
scenario.  The  reason  for  running  them  separately  at  first  was  to  see  how  each 
axis  compared  to  the  more  realistic  simulation  where  all  three  were  off  from  the 
orbital  frame.  Notice  on  all  simulations,  the  control  deflections  are  relatively  small, 
Sfi  <5°.  It  is  important  that  these  angles  remain  small  since  this  is  a  linearized 
model.  At  higher  altitudes  the  maximum  control  angle  deflections  actually  decrease 
since  the  controller  is  a  proportional  controller.  At  higher  altitudes  with  the  same 
initial  conditions,  the  rotational  rate  decreases  due  to  a  smaller  force  imparted  on 
the  control  panels. 


Control  Angle  Deflections 


(a) 


(b) 


Figure  5.6:  (a)  Angular  displacement  with  roll,  pitch  and  yaw  offset  10°. 

(b)  Control  angles  with  roll,  pitch  and  yaw  offset  10°. 


5.1.2  Altitudes  Above  200km.  Since  we  are  interested  in  determining  how 
altitude  affects  the  settling  time,  the  same  procedure  was  done  for  altitudes  of  300km, 
400km,  500km  and  600km.  See  Appendix  E  for  the  plots  not  shown  in  this  chapter. 
For  each  simulation,  the  time  constants  were  found  and  are  summarized  in  Tables 

5.3  and  5.4.  Table  5.3  is  the  case  where  each  axis  was  deflected  separately  and  Table 

5.4  is  the  case  where  all  three  axes  were  deflected  simultaneously.  In  either  case, 
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there  isn’t  much  difference  between  both  tables.  By  using  the  same  exponential 
curve  fit  algorithm  as  done  in  section  5. 1.1. 2,  an  exponential  curve  was  fit  to  the 
time  constants  in  Table  5.4  for  the  roll,  pitch  and  yaw  axes.  Figures  5.7,  5.8  and  5.9 
show  how  the  time  constants  grow  exponentially  with  altitude. 


Table  5.3:  Time  constants  for  each  axis  at  the  various  altitudes  where  roll,  pitch 
and  yaw  were  offset  separately. 


Altitude  (km) 

Time  Constant  (s) 

Roll  Axis 

Pitch  Axis 

Yaw  Axis 

200 

500 

130 

85 

300 

975 

1550 

950 

400 

6600 

10200 

6200 

500 

40000 

55000 

30000 

600 

275000 

265000 

260000 

Table  5.4:  Time  constants  for  each  axis  at  the  various  altitudes  where  roll,  pitch 
and  yaw  were  deflected  simultaneously. 


Altitude  (km) 

Time  Constant  (s) 

Roll  Axis 

Pitch  Axis 

Yaw  Axis 

200 

481 

130 

85 

300 

960 

1550 

950 

400 

6800 

10200 

6000 

500 

40000 

55000 

30000 

600 

265000 

265000 

260000 
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Figure  5.7:  Time  constants  for  the  roll  axis  increase  exponen¬ 
tially  with  increasing  altitude. 


Figure  5.8:  Time  constants  for  the  pitch  axis  increase  expo¬ 

nentially  with  increasing  altitude. 
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Figure  5.9:  Time  constants  for  the  yaw  axis  increase  expo¬ 

nentially  with  increasing  altitude. 


5.2  Pitch  as  the  Intermediate  Moment  of  Inertia  Axis 


The  satellite  would  be  in  its  most  stable  configuration  with  the  pitch  axis  as 
the  major  or  minor  MOI  axis.  That  may  not  always  be  the  case,  so  the  pitch  axis  was 
set  as  the  intermediate  axis  which  would  make  the  attitude  unstable  (see  equation 
5.3  for  the  MOI  matrix).  For  this  simulation,  theta  was  offset  by  1  degree  for  each 
axis,  the  altitude  was  set  to  400km  and  density  was  set  to  zero.  This  shows  the 
satellite’s  attitude  is  unstable  with  pitch  as  the  intermediate  axis  as  can  be  seen  in 
Figure  5.10. 


MOI  = 


A 

0 

0 

107 

0 

0 

0 

D 

0 

= 

0 

92 

0 

0 

0 

c 

0 

0 

68 

(5.3) 


When  the  atmospheric  density  is  put  back  into  the  simulation,  the  controller  is  able 
to  return  the  satellite  back  to  flying  at  its  intended  attitude.  See  figure  5.11 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  5.10:  Pitch  as  the  intermediate  MOI  axis  and  atmo¬ 
spheric  density  p  =  0. 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  5.11:  Pitch  as  the  intermediate  MOI  axis  and  atmo¬ 
spheric  density  shown  in  Table  5.1. 
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5.3  Roll,  Pitch,  and  Yaw  Having  the  Same  MOI 


Since  the  pitch  and  yaw  axes  have  the  same  proportional  constant  for  the 
controller  and  are  identical  in  dimensions,  they  should  have  the  same  time  constant. 
This  is  not  the  case  as  can  be  seen  in  Figure  5.12,  where  the  amplitude  of  the  yaw 
oscillations  end  up  being  smaller  than  the  pitch  axis.  This  is  due  to  the  roll  and  yaw 
axis  being  coupled.  The  yaw  and  roll  axes  work  work  together  thus  having  more 
control  authority.  Therefore  the  yaw  axis  decreases  in  amplitude  at  a  faster  rate 
than  the  pitch  axis. 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  5.12:  Angular  position  and  all  three  axes  having  the 

same  MOI. 


VI.  Conclusions  and  Recommendations 


6. 1  Conclusions 

To  be  able  to  launch  a  satellite  to  orbit,  and  have  the  spacecraft  orient  itself 
automatically  with  a  very  simple  attitude  control  system  is  very  desirable  and  has 
many  advantages.  A  simple  control  system  would  be  less  costly  and  less  prone 
to  failure  than  complex  systems  such  as  control  moment  gyros,  reaction  wheels  or 
thrusters.  Since  most  satellites  have  Earth  sensing  missions,  they  are  typically  nadir 
pointing  and  require  some  form  of  attitude  control.  These  satellites  typically  have 
their  body  reference  frame  aligned  with  the  orbit  reference  frame  and  are  usually  in 
LEO.  Satellites  are  put  in  LEO  for  various  reasons  such  as  payload  limitations  and 
cost  issues  and  have  to  compensate  for  atmospheric  effects.  Using  the  atmosphere  to 
control  a  satellite’s  attitude  was  investigated  and  the  results  presented  here  show  that 
it  is  a  feasible  alternative  to  the  more  expensive  and  complex  methods  mentioned 
earlier  for  satellites  in  LEO. 

The  atmospheric  density  decreases  exponentially  with  altitude  and  therefore 
an  atmospheric  drag  type  control  system  has  an  altitude  limitation  where  it  no 
longer  has  enough  control  authority  to  overcome  disturbance  torques.  The  effects  of 
increasing  altitude  can  be  seen  by  comparing  the  settling  time  at  different  altitudes. 
When  a  satellite’s  attitude  is  perturbed  from  its  alignment  with  the  orbital  reference 
frame,  it  takes  time  for  the  attitude  control  system  to  reduce  the  amplitude  of 
oscillations  and  bring  it  back  to  its  proper  attitude.  At  altitudes  of  200  km,  the 
settling  time  is  on  the  order  of  minutes  as  compared  to  a  600km  altitude  where 
the  settling  time  is  on  the  order  of  weeks.  Settling  times  are  also  affected  by  the 
satellites  physical  dimensions  and  mass.  An  increase  in  drag  panel  area,  lengthening 
the  control  arms  or  decreasing  mass  will  reduce  settling  time,  however,  the  settling 
times  would  still  be  approximately  the  same  order  of  magnitude.  The  linear  analysis 
performed  here  has  a  small  control  angle  limitation.  A  nonlinear  analysis  would 
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allow  the  control  panels  to  be  fully  extended  so  the  projected  areas  were  maximized 
which  will  also  reduce  settling  time. 

According  to  the  simulation  results,  the  oscillations  decay  exponentially  but 
will  never  quite  decay  to  zero.  In  reality,  there  is  a  limit  to  the  accuracy  of  the 
attitude  sensors  and  how  finely  the  drag  panel  angles  can  be  controlled,  therefore,  it 
can  be  expected  that  the  satellite  will  continue  to  oscillate  about  the  orbital  reference 
frame  at  plus  or  minus  some  small  angle.  Depending  on  the  pointing  requirements 
of  the  payload,  these  oscillations  may  not  be  acceptable.  If  so,  a  stabilized  payload 
could  be  added  to  the  spacecraft  to  counter  these  oscillations. 

The  bottom  line  is  that  atmospheric  drag  can  effectively  be  used  to  control  a 
spacecraft  attitude  in  LEO.  Disturbances  in  attitude  can  be  overcome  resulting  in  a 
stable  attitude.  The  spacecraft’s  control  system  would  consist  of  a  less  costly,  simple 
control  system  with  simple  actuators  less  prone  to  failure.  The  research  presented 
here  is  the  foundation  for  future  study  which  can  eventually  lead  to  an  actual  flight 
test  of  a  drag  controlled  satellite. 

6.2  Recommendations  for  Future  Research 

The  research  presented  here  can  be  very  useful  in  several  areas  of  continued 
research.  Future  research  should  be  built  upon  eventually  leading  to  the  flight  test 
of  an  actual  aerodynamically  controlled  satellite.  Below  are  several  ideas  to  build 
upon  this  research. 

6.2.1  Higher  Fidelity  Nonlinear  Computer  Model.  The  next  step  with 
this  topic  would  be  to  add  estimation  algorithms  into  the  model  and  eventually 
do  a  nonlinear  analysis.  A  nonlinear  computer  model  would  better  represent  the 
dynamics  of  an  actual  aerodynamically  controlled  satellite  system.  The  final  step 
would  be  to  set  the  spacecraft  in  a  tumble  with  the  nonlinear  model  and  see  if 
it  can  recover.  Other  factors  should  be  added  to  the  model  such  as  a  rotating 
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atmosphere,  elliptical  orbits,  non-equatorial  orbits,  drag  panels  with  mass,  random 
external  torques  from  solar  winds,  oblate  Earth,  and  an  orbit  decay  algorithm  to 
determine  the  approximate  orbit  lifetime  for  a  satellite  of  a  given  mass,  projected 
area  and  altitude.  This  way,  the  optimum  orbit  could  be  selected  based  on  lifetime 
and  payload  capabilities/requirements.  From  this  point,  satellite  dimensions  could  be 
determined  for  the  most  effective  control.  The  altitudes  where  the  control  authority 
is  not  enough  to  overcome  perturbations  can  also  be  determined.  Note,  orbit  decay 
calculations  are  only  approximate  since  there  are  so  many  random  variables  that 
affect  the  atmospheric  density. 

6.2.2  Solar  Drag  Panels.  Since  satellites  usually  have  solar  panels  to 
generate  electrical  power,  the  drag  panels  should  be  replaced  with  solar  panels.  In 
this  research,  the  mass  of  the  panels  was  neglected  however,  a  solar  panel  mass  should 
be  included  in  the  model.  Including  the  mass  of  the  panels  in  the  model  would  also 
require  some  vibrational  analysis  to  determine  potential  problems  with  oscillations 
of  the  panels.  As  the  panels  are  being  moved,  toques  will  have  an  effect  on  the 
satellites  attitude  and  should  be  accounted  for.  Finally,  research  the  feasibility  of 
tracking  the  sun  with  one  or  more  panels  while  maintaining  attitude. 

6.2.3  Other  Related  Research.  Some  students  have  done  research  on  main¬ 
taining  satellite  formations.  This  model  could  be  used  with  their  research  to  control 
both  the  satellites  attitude  and  station  keeping  at  the  same  time.  Other  designs 
using  atmospheric  drag  for  attitude  control  should  be  considered  and  compared  with 
this  and  other  research  to  determine  the  best  and  most  cost  effective  design.  Finally, 
research  should  be  done  to  determine  if  drag  panels  could  be  used  for  momentum 
dumping  of  reaction  wheels. 
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Appendix  A.  MathCad  Worksheets 


A.l  MathCad  Nonlinear  Torque  Calculation  Worksheet 
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First  find  Rba  using  the  definition  and  by  using  3-2-1  Euler  Angle  rotations 
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\0  — sin(O)  cos(0) ) 
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cos(O)  sin(O)  0  10  -sin(0)  cos(0) 
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Calc  projected  areas  in  the  A  frame: 
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Calc  Forces  in  the  A  frame  for  the  Ailerons: 


£ 


X 

00 

00 

Os 


II 

Mh 

p 


X 

00 

00 

Os 


<5 


79 


rightF  :=  (°  1  °)  %A.right  FnghtA  (°  1  °>  FnghtF  =  -2.988  x  10 


=  idoxnviJ  ( i  0  0).vdoinviJ-d°invw3i  v%.(  x  0  0)  =:  id°inv7 
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Convert  Ailerons  to  B  frame  from  F: 
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r  is  the  position  vector  from  the  s/c  eg  to  the  centroid  of  the  panels  B  Frame: 
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Appendix  B.  Matlab®  Code 

B.l  Matlab®  Eigenvalues  and  Eigenvectors  with  Torques 

Listing  B.l:  This  Matlab®  code  calculates  the  eigenvalues  and  eigenvectors  for 

the  model  validation. 

(appendix2/EigsEOMw  Torques. m) 

7.  D  .  Guettler 
clc 

clear  all 
Alt  =300;  7.  km 

5  syms  A  B  C  n  Rho  Vel 

7.  MAT1  =  [0  00100;  000010;  000001;... 

7,  (  (C-B) /A)  *n~2  0  0  0  0  (  (C  +  A-B)  /A)  *n  ;  .  .  . 

7.  0  -4.621*Rho*Vel~2/B  0  0  0  0;... 

10  7.  0  0  ((A-B)*n~2-4.414*Rho*Vel~2)/C  ( (B-A-C) /C)  *n  0  0] 

7.MAT2  is  the  4x4  matrix  of  MAT1 
MAT2  = [0  010;  0001;... 

( (C-B) /A) *n~2  0  0  ( ( C  +  A-B ) / A ) *n  ;  .  .  . 

15  0  ((A-B)*n~2-4.414*Rho*Vel~2)/C  (  (B-A-C) /C)  *n  0] 

7.  Orbital  rate  for  300  km  altitude  in  rad /sec 

[n  ,  T  ,  f  req  ,  Vel  ]  =  OrbRate(Alt)  7.n  =  0.001157  7«for  300km  alt 

20  7. Calculate  Density 
Rho  =  atmos.exp ( Alt ) 

7. Rho  =  0  ; 

7.M0I  B  is  always  major  axis  ,  C  is  always  minor  axis 
25  A  =  20  . 

B  =  30  . 

C  =  10  . 


7.  MAT1  =  [0  00100;  000010;  000001;... 

30  7.  (  (C-B) /A) *n~2  0  0  0  0  (  (C  +  A-B)  /A) *n ; . . . 

7.  0  -4.621*Rho*Vel~2/B  0  0  0  0;... 

7.  0  0  (  (A-B)  *n  ~  2 -4 .414*Rho*Vel~2)  /  C  (  (B-A-C) /C)  *n  0  0] 


7.MAT2  is  the  4x4  matrix  of  MAT1 
35  MAT2  = [0  010;  0001;... 

( (C-B) /A) *n~2  0  0  ( ( C  +  A-B ) / A ) *n  ;  .  .  . 

0  ((A-B)*n~2-4.414*Rho*Vel~2)/C  ( (B-A-C) /C) *n  0] 


7.  7. From  6x6  Matrix 
40  7.  [EigVecl  ,  EigVal  1  ]  =  e ig  ( MAT1  )  ; 
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l  EigVall  =  simplify (EigVall ) ; 
V.  EigVall  =  vpa(EigVall  ,3) 

1  EigVecl  =  simplify (EigVecl )  ; 
7.  EigVecl  =  vpa(EigVecl  ,3) 

'/.From  4x4  matrix 
[EigVec2 , EigVal2] =eig (MAT2 ) ; 
EigVal2  =  s impl if y ( EigVal2 ) ; 
EigVal2  =  vpa ( EigVal2 , 3) 
EigVec2  =  s impl if y ( EigVe c2 ) ; 
EigVec2  =  vpa ( EigVec2  ,  3) 


'/.Find  ICs  for  each  run 
EigValCasel=EigVal2 (1,1) 

ICl  =  (l+i) * [EigVec2 (1  ,  1)  EigVec2(2,l)  EigVec2(3,l)  EigVec2 (4 , 1)  . .  . 
]’  +  ... 

( 1  —  i ) * [EigVec2 (1 ,2)  EigVec2(2,2)  EigVec2(3,2)  EigVec2 (4 , 2) ]  ’ 


EigValCase2=EigVal2 (3,3) 

IC2  =  (l+i) * [EigVec2 (1 ,3)  EigVec2(2,3)  EigVec2(3,3)  EigVec2 (4 , 3)  . .  . 
]’  +  ... 

( 1 - i ) * [ EigVec2 ( 1 , 4)  EigVec2(2,4)  EigVec2(3,4)  EigVec2 (4 , 4) ]  ’ 


B.2  Matlab®  Curve  Fit  Algorithm 

This  algorithm  uses  a  peak  finding  function  to  locate  all  the  peaks  of  the  os¬ 
cillations  of  the  angles  between  the  orbit  and  body  frame.  Once  the  peaks  were 
located,  the  main  code  fit  an  exponential  curve  to  the  points  based  on  an  approxi¬ 
mated  time  constant.  The  time  constant  was  then  varied  until  an  exponential  curve 
closely  matched  the  peaks.  This  code  was  used  to  fit  an  exponential  curve  to  the 
decaying  oscillations  for  the  time  constant  calculations  in  Chapter  V. 

Listing  B.2:  This  Matlab®  code  fits  the  peak  amplitudes  to  an  exponentially 

decaying  function. 

(appendix2  /  TimeConst  Calc .  m) 

7,  peaks  =  function  peakf  inder  (  roll  ,  t  ime  ) 

7,  D  .  Guettler 

clear ; clc  ; 

load  thetaout . mat 

7.Pick  tau  as  time  constant  compare  plot 
tau=38000 
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“i  ROLL 

70  peaks  =  peakfinder(theta(2,:)’,theta(l,:)’); 
l  PITCH 

70  peaks  =  peakfinder(theta(3,  :)  ’  ,  thet  a  ( 1  ,  :)  ’  )  ; 

°i  YAW 

peaks  =  peakfinder(theta(4,:)’,theta(l,:)’); 

7» Create  vectors  for  time  calculation 

t=peaks ( : , 1) 
y=peaks ( : , 2) 

70  ere  ate  plot 

X  =  [ exp ( -t / tau ) ] ; 
i  Calculate  model  coefficients 
a  =  X\y 

T  =  (0:1: 200000)  ’ ; 

Y  =  [exp (-T/tau) ] *a ; 

i  Create  figure 
figurel  =  figure; 
i  Create  axes 

axesl  =  axes (’ Parent ’, f igure 1 ) ; 
box ( 1  on  ’  )  ; 
grid ( ’ on  ’ )  ; 
hold ( ’ all  ’  )  ; 
i  Create  plot 

plotl  =  plot(T,Y, ’Parent’, axesl),  grid  on 
i  Create  plot 

plot2  =  plot (t , y ,  ’LineStyle’  ,  ’none1  ,  ’Marker1  ,  ’o’  ,  .  .  . 

1  Parent ’  , axe  s 1 )  ; 

70  Create  x  lab  el 
xlabel(’Time  (s)1); 

°i  Create  y  lab  el 

ylabel (’ Amplitude  (rad)’); 

I  Create  legend 

legendl  =  legend( axesl, {’Exponential  Curve  Fit’,... 

’Peak  Locations ’},  ’Position’  ,  .  .  . 

[0.5827  0.8245  0.3232  0.1002]); 

Listing  B.3;  This  Matlab®  code  determines  the  positive  peak  amplitudes  of  oscil¬ 
lation. 

(appendix2 /PeakF  inder .  m) 
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“/»  peakf  inder  .  m 
%  D .  Guettler 

function  peaks  =  peakf inder ( roll , t ime ) 

'/.Need  roll,  time  as  column  vectors  of  equal  dimension 
"/.Includes  initial  and  final  points 

[cols ,rows]=size(roll) ; 


if  roll ( 1 , 1 ) >0 

r ollpeaknums ( 1 ) =1 ; 

j=2; 

else  j  =1 ; 
end 

for  i =2 : ( cols -1 ) 

if  roll(i,l)>roll(i-l,l)  &  roll(i,l)>roll(i+l,l)... 
&  roll ( i  ,  1 )  >0 
r ollpeaknums ( j ) =i ; 

j=j+i ; 

end 

end 

r ollpeaknums (j)=cols; 
for  i =1 : j 

peaks (i , 1) =time ( r ollpeaknums (i) ,1) ; 
peaks (i ,2) =roll ( r ollpeaknums ( i )  ,  1 )  ; 

end 
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Appendix  C.  The  Atmosphere 

C.l  Altitude  vs  Density  Plot  (NASA)  [3,  page 

i - 1 - r 


1000 


900 


800 


700 


-  600 
E 


500 


400 


300 


200  - 


100 


i - 1 - 1 - 1 - r 


^ — i — i — r 


Upper  Limit  of  ^ 
COESA  (1958) 
Standard 


COESA  (1976) 
COESA  (1962) 


- Minzner  at  at.  (1959) 

(COESA  (1958) 

(Minzner  and  Ripley  (1956) 
ooooooo  Rocket  Panel  (1952) 
o-o-c-o  Warfield  (1947) 


10' 1 


J _ I _ I _ L 

10'1 


10"  0  10' 


Atmospheric  Density  (kg  ni  ) 

Figure  C.l:  Atmospheric  Density/ Altitude  plot 


C.2  Matlab®  Density  Calculations) 

For  the  computer  simulations  (see  Appendix  C),  an  exponential  atmospheric 
model  was  used.  This  model  was  tested  at  the  following  altitudes  and  compared 
to  the  density  plot,  Figure  C.l.  The  results  are  displayed  in  Table  C.l  exponential 
atmospheric  density  calculations  in  Matlab®  compare  with  those  found  in  the  plot. 


Table  C.l:  Altitude  and  Density  Results  from  Matlab®  . 


Altitude 

Calculated  Density 

Density  from  Figure  C.l 

150 

2.07xl0~9 

2xl0~9 

200 

2.79xl0_1° 

2xlO“10 

250 

7.25xl0-11 

7xl0”n 

300 

2.42xl0”n 

lxl0~n 

350 

9.52xl0~12 

7xl0~12 

400 

3.73xl0~12 

3xl0“12 

450 

1.59xl0-12 

lxl0~12 

500 

6.97xl0-13 

5xl0~13 

550 

3.18xl0-13 

2xl0~13 

600 

1.45xl0“13 

lxlO”13 
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Appendix  D.  Simulink®  Model 

D.l  Simulink®  Model  Top  Level 

The  following  figures  illustrate  the  computer  model  used  to  simulate  the  dy¬ 
namics  of  using  atmospheric  drag  for  attitude  control.  The  software  used  is  Simulink® 
Figure  D.l  is  the  overall  Simulink®  model.  It  is  laid  out  such  that  the  inputs  are 
along  the  left  side,  the  outputs  are  on  the  right.  This  model  represents  equation 
3.47  with  all  the  variables  on  the  input  side.  The  model  outputs  the  control  angles, 
angular  accelerations,  angular  rates  and  attitude  angle  with  respect  to  the  orbital 
frame.  All  satellite  dimensions  are  input  as  meters,  angles  are  input  as  degrees  and 
altitude  is  input  as  km.  All  the  outputs  containing  angles  are  in  degrees. 
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Figure  D.l:  Overall  Simulink®  Model 
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D.2  Function:  Orbit,  m 


Listing  D.l:  This  Mat  lab®  function  calculates  the  approximate  density  of  the 

atmosphere  at  a  given  altitude,  tangential  orbital  velocity  for  a  circular  orbit,  orbital 
rate,  orbital  period  and  orbital  frequency. 

(appendix4/Orbit.m) 


function  [rho,  Vel , n , T , freq] = Orbit ( Alt ) 

%  Exponential  Atmospheric  Model 
7o  Accepts  as  INPUT  the  orbital  radius  of  the 
7  satellite;  OUTPUTS  the  atmospheric  density 
7  for  that  altitude  based  on  the  mean  equatorial 
7  radius  of  Earth  (density  units:  kg/m~3). 

7  NOTE:  really  only  valid  for  altitudes  less  than 
7  1000  km.  Taken  from  Vallado  ,  pg  537. 

7.  Distances  in  km 


7  Calculate  Orbit  properites  (D .  Guettler) 

mu  =  3.986*10~14  ;  7 ( km "3/ s ~2)  Standard  Gravitational 

'/.Parameter 

REarth  =  6378135  ;  7(km)  Radius  of  Earth  in 

ROrbit  =  REarth  +  Alt  *  1000  ;  ‘/.Radius  from  Earth’s 

%  center  to  orbit 


OrbRate  =  sqrt (mu/ROrbit ~3)  ; 

Vel  =  sqrt (mu/ROrbit )  ; 
n  =  OrbRate  ; 
T  =  2* pi  * ( 1 / OrbRat e )  ; 
freq=l/T  ; 


%  OrbRate  is  the  orbital 
'/.rate  in  rad/sec 
7,  Orbit  Velocity 
/{Orbital  rate  (rad/sec) 
7.T=orbital  (s) 

7»  Frequency  (orbits/ sec) 

7.  Calculate  Density  (B.  Hajovsky) 

7.  Term  Definitions: 

7.  hO  =  Base  Altitude  (km) 

7«  rhoO  =  Nominal  Density  (kg/m~3) 

7«  H  =  Scale  Height  (km) 

7«  Define  the  Earth  radius  to  calculate  the  altitude: 
R=Alt ; 
r  =  R ; 
rhoO  =0 ; 
hO  =0  ; 

H  =  0  ; 

if  r  >  =  0  &  r  <25 


hO  =  0; 
rhoO  =  1.225; 
H  =  7.249; 
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elseif  r>=25  &  r<30 
hO  =  25; 

rhoO  =  3 . 899e-2 ; 

H  =  6 . 349 ; 
elseif  r>=30  &  r<40 
hO  =  30; 

rhoO  =  1 . 774e-2 ; 

H  =  6 . 682 ; 
elseif  r>=40  &  r<50 
hO  =  40; 

rhoO  =  3.972e-3; 

H  =  7.554; 
elseif  r>=50  &  r<60 
hO  =  50; 

rhoO  =  1 . 057  e  -3  ; 

H  =  8.382; 
elseif  r>=60  &  r<70 
hO  =  60; 

rhoO  =  3 . 20  6  e  -4  ; 

H  =  7.714; 
elseif  r>=70  &  r<80 
hO  =  70; 

rhoO  =  8 . 770  e -5 ; 

H  =  6 . 549 ; 
elseif  r>=80  &  r<90 
hO  =  80; 

rhoO  =  1 . 905e-5 ; 

H  =  5 . 799 ; 

elseif  r>=90  &  r<100 
hO  =  90; 

rhoO  =  3 . 396e-6 ; 

H  =  5 . 382  ; 

elseif  r>=100  &  r<110 
hO  =  100; 
rhoO  =  5 . 297  e  -7  ; 

H  =  5 . 877 ; 

elseif  r>=110  &  r<120 
hO  =  110; 
rhoO  =  9.661e-8; 

H  =  7.263; 

elseif  r>=120  &  r<130 
hO  =  120; 
rhoO  =  2 . 438  e  -8  ; 

H  =  9.473; 

elseif  r>=130  &  r<140 
hO  =  130; 
rhoO  =  8 . 484e -9 ; 

H  =  12.636; 
elseif  r>=140  &  r<150 
hO  =  140; 
rhoO  =  3 . 845  e -9 ; 
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H  =  16 . 149 ; 
elseif  r>=150  &  r<180 
hO  =  150; 
rhoO  =  2.070e-9; 

H  =  22.523; 
elseif  r>=180  &  r<200 
hO  =  180; 
rhoO  =  5.464e-10; 

H  =  29.740; 
elseif  r>=200  &  r<250 
hO  =  200; 
rhoO  =  2.789e-10; 

H  =  37.105; 
elseif  r>=250  &  r<300 
hO  =  250; 
rhoO  =  7 . 248  e  -11; 

H  =  45.546; 
elseif  r>=300  &  r<350 
hO  =  300; 
rhoO  =  2.418e-ll; 

H  =  53.628; 
elseif  r>=350  &  r<400 
hO  =  350; 
rhoO  =  9.518e-12; 

H  =  53.298; 
elseif  r>=400  &  r<450 
hO  =  400; 
rhoO  =  3.725e-12; 

H  =  58 . 515  ; 
elseif  r>=450  &  r<500 
hO  =  450; 
rhoO  =  1 . 585  e  -  1 2  ; 

H  =  60 . 828; 
elseif  r>=500  &  r<600 
hO  =  500; 
rhoO  =  6.967e-13; 

H  =  63.822; 
elseif  r>=600  &  r<700 
hO  =  600; 
rhoO  =  1.454e-13; 

H  =  71 . 835  ; 
elseif  r>=700  &  r<800 
hO  =  700; 
rhoO  =  3.614e-14; 

H  =  88.667; 
elseif  r>=800  &  r<900 
hO  =  800; 
rhoO  =  1.170e-14; 

H  =  124.64; 

elseif  r>=900  &  r<1000 
hO  =  900; 
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rhoO  =  5 . 245  e  -  1 5  ; 

H  =  181.05; 
elseif  r >1000 
hO  =  1000; 

150  rhoO  =  3.019e-15; 

H  =  268 . 00  ; 

else 

end 

155  rho  =  rhoO *  exp (- (r-hO ) /H )  ; 
"/,  rho=0  ; 


D.3  Function:  Constants,  m 

Listing  D.2:  This  Matlab®  function  calculates  the  constants  in  equations  3.46  and 
also  calculates  the  MOI  matrix  based  on  a  given  satellite  mass. 
(appendix4/Constants.m) 

function  [CTArm  ,  CBArm  ,  CLArm  ,  CRArm  ,  CLAil  ,  CRAil  ,  InvMOI  ,  .  .  . 

CTheta2 , CTheta3 , MOI ]  =  Constants  (SatDims) 

"/.  This  block  calculates  the  constants  for  the 
"/.  linearized  torques  and  contains  all  the  satellite 
5  "/.  dimensions.  See  MathCad  File  for  explanation. 

"/.  D  .  Guettler 

"/,  Del  care  Variables 

Lpanel=  SatDims  ( 1)  ;  "/.Length  of  drag  panel 

'/.(front  to  back  when  retracted)  (m) 

Wpanel=  SatDims  (2)  ;  "/.Width  of  panel  (m) 

LArm  =  SatDims  (3)  ;  "/.Length  of  control  arm  (m) 

15  PhiOArm  =  SatDims  (4)  *pi  / 180  ;  "/.Linearization  angle  for 

"/.arms  (radians) 

LAil  =  SatDims  (5)  ;  "/.Length  of  Aileron  spanwise  (m) 

WAil  =  SatDims  (6)  ;  "/.Width  of  Aileron  chordwise  (m) 

Lbl  =  SatDims  (7)  ;"/,  cubesat  length  along  bl  axis  (m) 

20  Lb2  =  SatDims  (8)  ;  “/.Cubesat  width  along  b2  axis  (m) 

Lb3  =  SatDims  (9)  ;  "/.cubesat  Height  along  b3  axis  (m) 

PhiWedge  =  SatDims  ( 10)  *pi /180 ;  "/.one  half  the  total 

"/.wedge  angle  (Radians) 

Mass  =  SatDims  (11);  "/.Mass  of  sat  in  kg 
25  "/.(used  to  calc  MOI  matrix) 

CD  =  SatDims  (12);  "/.Drag  Coefficient 
Rho  =  SatDims (13); 

V  =  SatDims ( 14) ; 

"/.  Calculate  Constants 
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•/.  Const  mult  iplying  deltaPhiTopArm 

CTArm=(.5*  Lpanel “ 2* Wpanel *  cos ( PhiOArm ) * s in ( PhiOArm )  ... 

35  +  Lpanel * Wpanel * LArm* cos ( PhiOArm )* sin ( PhiOArm )  ... 

+  . 75* Lpanel * Wpanel * Lbl  *  cos (PhiOArm) ~2* sin (PhiOArm)  ... 

-  . 25* Lpanel * Wpanel * Lbl * s in ( PhiOArm )  ... 

+  . 75* Lpanel * Wpanel * Lb3 *  cos ( PhiOArm )  ... 

-  ,75*Lpanel*Wpanel*Lb3*cos( PhiOArm ) ~3) *CD*Rho*V~2; 

40  /  Const  multiplying  deltaPhiBot Arm 

CBArm=-CTArm ; 

•/.  Const  mult  iplying  deltaPhiLeftArm 

CLArm  =  (.25  * Lpanel * Wpanel *Lb 1  *  sin ( PhiOArm )  ... 

-. 75*Lpanel*Wpanel*Lb2*cos (PhiOArm) . . . 

45  +. 75*Lpanel*Wp anel* Lb 2*cos (PhiOArm) "3  .  .  . 

- 1  * Lpanel * Wpanel * LArm* cos ( PhiOArm )* sin ( PhiOArm )  .  .  . 

-.  5* Lpanel ~2* Wpanel *  cos ( PhiOArm )* sin ( PhiOArm )  .  .  . 

-.75*Lpanel*Wpanel*Lbl*cos (PhiOArm) ~2*  s in ( PhiOArm ))*CD*Rho*V... 
-2; 

7,  Const  mult  iplying  deltaPhiRightArm 
50  CRArm  =  -CLArm  ; 

7,  Const  mult  iplying  deltaPhiLeft  Aileron 

CLAi 1  =  ( 1 . 5* LAil ~ 2* WAil *  cos ( PhiWedge ) " 2  * s in ( PhiWedge )  ... 

- . 5* LAil * WAi 1  * Lb2 * s in ( PhiWedge )  .  .  . 

+ 1 . 5*  LAil  *  WAil  *Lb2  *  cos  (  PhiWedge  )  ”  2  *  s  in  (  PhiWedge  )  ... 

55  -,5*LAil~2*WAil*sin(PhiWedge))*CD*Rho*V~2; 

7.  Const  mult  iplying  de  It  aPhiRight  Aileron 
CRAil=-CLAil ; 

7.  Const  mult  iplying  Theta  2  term 

CTheta2  =  (-l  * Lpanel *Wpanel  *  Lb  1  *  cos ( PhiOArm ) “ 2* s in ( PhiOArm )  ... 

60  -1  *Lpanel ~2*Wpanel*cos (PhiOArm) *sin (PhiOArm)  ... 

+  1  * Lpanel * Wpanel *Lb3 *  cos ( PhiOArm ) ~3  ... 

-2  * WAi 1 “ 2* LAil *  cos ( PhiWedge )* sin ( PhiWedge )  ... 

-2  * Lpanel *Wpanel * LArm  *  cos ( PhiOArm )* s in ( PhiOArm )  ... 

+  2  * WAi 1 “ 2* LAil *  cos ( PhiWedge )“ 3* s in ( PhiWedge )  ... 

65  -1  * Lpanel * Wpanel *Lb3 *  cos ( PhiOArm ))* CD *Rho *V “ 2 ; 

7.  Const  mult  iplying  Theta  2  term 

CTheta3  =  (-l  * Lpanel " 2* Wpanel *  cos ( Phi OArm )* s in ( Phi OArm )  ... 

-2  * Lpanel * Wpanel * LArm *  cos ( PhiOArm )* sin ( PhiOArm )  ... 

+  1  * Lpanel * Wpanel *Lb2 *  cos ( PhiOArm ) ~3  ... 

70  -1  * Lpanel * Wpanel *Lb 1  *  cos ( PhiOArm ) ~ 2* s in ( PhiOArm )  ... 

-1  *Lpanel*Wpanel*Lb2*cos (PhiOArm) ) *CD*Rho*V~2; 

7.  Calc  MOI  Where  pitch  is  major,  Roll  is  intermediate, 

75  7,  Yaw  is  minor 

Length  =  1 ; 

Width  =  . 8 ; 

Height  =  1.25; 

80  Ir oil  =  ( 1 / 1 2 ) *Mass * ( Width ~ 2  +  He ight  “  2)  ; 
Ipitch=(l/12)*Mass*(Length"2+Height~2) ; 
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Iy aw =(1/12) *Mass*(Length~2+Width"2) ; 

7oXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

"/0MOI=[Iroll  0  0;0  Ip itch  0 ;  0  0  Iyaw];  7«(kg*m~2) 
MO  I  =  [10  0  0  ;  0  20  0  ;  0  0  30]  ; 

“/oXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

7°  Calc  Inv  MO  I 
InvMOI =MOI “ - 1 ; 


D.4  Function:  ControlTorques.m 

Listing  D.3:  This  Matlab®  function  calculates  Cx 3  (equation  3.50)and  multiplies 
it  with  the  deflections  of  the  control  surfaces. 

(appendix4/ControlTorques.m) 

function  [Angles , CX3dPhi ]  =  ControlTorques (CTArm , CBArm , CLArm , . . . 

CRArm , CL Ail , CRAil , InvMOI , C , ThetaDot , Theta) 

7«  This  block  calculates  the  torques  a  satellite 
7«  experiences  when  in  LEO,  from  Theta  (the  angles 
7.  between  the  B-frame  and  the  A-frame)  ,  PhiArms 
7.  (the  angle  of  each  control  arm  from  the  sat 
I  body)  and  PhiFlaps  (the  angle  of  the  flaps 
7.  rotating  about  the  arms)  .  it  is  assumed  the 
7.  inward  pointing  normal  direction  for  the  panels 
7.  never  sees  the  wind. 

%D .  Guettler 

'/.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
7«EIG  TEST  Comment  out  for  control 
7.  RollContr  ol  =0  ; 

7.  PitchControl=0; 

7.  YawControl  =0  ; 

7«  Uncomment  for  control 
RollControl=C(l)*ThetaDot (1) ; 

PitchControl=C(2)*ThetaDot (2) ; 

YawControl=C (3) *ThetaDot (3) ; 

'/.xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
dPhiRightAil=RollControl ; 
dPhiLef t Ail=-RollControl ; 

dPhiTopArm=-PitchControl ; 
dPhiBotArm=PitchControl ; 

dPhiLef t Arm=YawControl ; 
dPhiRight Arm=- YawControl ; 

Angles =[ dPhiRight Ail  dPhiLeftAil  dPhiTopArm  dPhiBotArm . . . 
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dPhiLeftArm  dPhiRightArm] *180/ pi ’ ; 

35  CX3dPhi=InvM0I*  [ CLAil * dPhiLef  tAil+CRAil*dPhiRightAil  ; 
CTArm*dPhiTopArm+CBArm*dPhiBotArm; 
CLArm*dPhiLeftArm+CRArm*  dPhiRightArm]  ; 


D.5  Function:  CalcMat.m 

Listing  D.4:  This  Matlab®  function  calculates  the  matrices  C\1  and  Cx1  in  equa¬ 
tions  3.48  and  3.49. 

(appendix4/CalcMat.m) 

function  [CX1,CX2]  =  Cal cMat ( MOI , n , CThet a2 , CThet a3 ) 

°/0  This  block  calculates  the  matrix  mult  ipl  i  ed  by 

°/0  the  ThetaDot  vector.  See  the  help  menu  for  details. 

"i  Calc  Inv  MOI 
InvMO I =M0I ~ -1 ; 

°/o  Calculate  the  matrix  that  multiplies  the  ThetaDot  Vector 

CXI  =  InvMO  I  *  [0  0  (MOI  (3 ,3)  +M0I  (1  ,  1) -MOI  (2 , 2)  )  ; 

000;  (MOI (2 ,2) -MOI (1 , 1) -MOI (3 ,3) )  0  0]*n; 

15  \ 

7,  Calculate  the  matrix  that  mult  iplies  the  ThetaDot  Vector 

CX2  =  InvMOI  *  [(MOI  (3 ,3)  -MOI  (2 , 2)  )  *n"2  0  0; 

20  0  CTheta2  0; 

0  0  (MOI (1 , 1) -MOI (2 ,2) )*n  ~  2  +  CThet  a3 ]  ; 


D.6  Function:  InitCond.m 

Listing  D.5:  This  Matlab®  function  inputs  the  initial  conditions  that  were  calcu¬ 
lated  from  the  eigenvector  analysis  for  model  validation.  Once  the  model  has  been 
validated,  these  can  be  turned  off  so  the  actual  simulations  can  be  run. 
(appendix4/InitCond.m) 

function  [ThetaDot , Theta]  =  InitCond(n) 

%  this  block  is  to  test  the  model  by  the  e igenval 
7„  eigenvec  method.  See  the  help  menu  for  details. 
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"/.  D  .  Guettler 

“/.  ******For  the  Linear  no  torques  case****** 

“/.  For  the  cases  where  A  =  20  B  =  10  C  =  30  using  the 
“/.  4x4  matrix  to  get  eigvecs 

'/.Case  1:  The  first  pair  of  complex  conj  eigvecs 

7.  Thet  a  =  [2  0  2]  ’  ; 

7.  ThetaDot  =  [2*n  0  -2*n]  ’  ; 

“/.Case  2  The  second  pair  of  complex  conj  eigvecs 
"/.Thet  a  =  [1 . 734  0  2]’; 

"/.  Thet  aDot  =  [n  0  -1.156*n]’; 

"/.  ******For  the  Linear  no  torques  case****** 

7«  For  the  cases  where  A  =  10  B  =  20  C  =  30  using  the 
7.  4x4  matrix  to  get  eigvecs 

'/.Case  1:  The  first  pair  of  complex  conj  eigvecs 

“/.Thet  a  =  [1  0  1]  *  ; 

'/.ThetaDot  =  [n  0  -n]  1  ; 

'/.Case  2  The  second  pair  of  complex  conj  eigvecs 
Theta= [-3  00]’; 

ThetaDot=[0  0  n] ’ ; 

"/.Thet  a  =  [2  0  0]  ’ ; 

“/.ThetaDot  =  [.  00232  0  0]  ’  ; 

“/0Theta=  [0  0  2]  »  ; 

“/.ThetaDot  =  [0  0  .00232]“; 

“/0Theta=  [1 .414  0  1 .414]  ’  ; 

“/.ThetaDot  =  [.  00163599  0  -.00163599]’; 

“/. Theta  =  [  -  1 .3101 .512]  ’  ; 

“/.ThetaDot  =  [.000874  0  .101]’; 

"/.  For  the  cases  where  A  =  20  B  =  10  C  =  30  using  the 
7.  6x6  matrix  to  get  eigvecs 

“/.Case  1:  The  first  pair  of  complex  conj  eigvecs 

7. Theta  =  [2  0  2]  ’  ; 

"/.ThetaDot  =  [2*n  0  -2*n]  ’  ; 

7. Case  2  The  second  pair  of  complex  conj  eigvecs 
7. Theta=[-3/n  0  3.46/n]’; 

“/.ThetaDot  =  [1.73  0  2]  ’  ; 

7. Theta  =  [0  2  2]  ’  ; 

7. Thet aDot  =  [ -2*n  0  0]’; 


"/.For  the  linear  with  torques  case 

"/.A  =  20  ,  B  =  30,  C  =  10  and  using  the  4x4  matrix  to  get  eigs 
“/.Case  1  (first  pair  of  complex  conj  eigvecs) 

7.  Thet  a  =  [2  0  0]  ’  ; 

“/.ThetaDot  =  [.  00232  0  0]’; 

“/.Case  2  (2nd  pair  of  complex  conj  eigvecs) 

7.  Thet  a  =  [0  0  2]  ’  ; 

“/.ThetaDot  =  [0  0  .0506]’; 
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Appendix  E.  Simulation  Plots 


E.l  Altitude  Variation  plots 

This  Appendix  contains  all  the  plots  for  angular  displacement,  angular  rate, 
angular  acceleration,  and  control  angles  for  the  simulations  determine  the  effects 
altitude  has  on  attitude  stabilization.  All  inputs  are  summed  up  in  the  following 
tables. 

Table  E.l:  Density,  velocity  and  orbital  rate  calculated  from  Simulink®. 


Altitude  (km) 

Calculated  Density  (®|) 

Velocity  (™) 

Orbital  Rate  (^) 

200 

2.79xlO~10 

7784 

0.001183 

300 

2.42xl0~n 

7726 

0.001157 

400 

3.73xl0~12 

7669 

0.001131 

500 

6.97xl0-13 

7613 

0.001107 

600 

1.45xl0~13 

7558 

0.001083 

Table  E.2:  Variables  input  into  the  Simulink®  model  for  all  the  altitude  variation 
simulations. 


Variable 

Value 

L  Panel 

1  m 

W Panel 

1  m 

L  Arm 

Irn 

fi^Arm 

45o 

Lau 

1  m 

wAil 

1  m 

1  m 

-^62 

1  m 

Tfcg 

1  m 

4*  W  edge 

22.5° 

Mass 

500% 

cD 

2.2 

A 

0 

0 

92 

0 

0 

0 

B 

0 

= 

0 

107 

0 

0 

0 

c 

0 

0 

68 

kg  ■  m2 


(E.l) 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 
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Figure  E.l:  Angular  displacement  at  200km  with  roll  offset 

10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.2:  Angular  rate  at  200km  with  roll  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.3:  Angular  Acceleration  at  200km  with  roll  offset 

10°. 
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Control  Angle  Deflections 
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Figure  E.4:  Control  angle  deflections  200  km  with  roll  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.5:  Angular  displacement  at  200km  with  pitch  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.6:  Angular  rate  at  200km  with  pitch  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.7:  Angular  Acceleration  at  200km  with  pitch  offset 
10°. 


Control  Angle  Deflections 


Figure  E.8:  Control  angle  deflections  200  km  with  pitch  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.9:  Angular  displacement  at  200km  with  yaw  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.10:  Angular  rate  at  200knr  with  yaw  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.ll:  Angular  Acceleration  at  200km  with  yaw  offset 
10°. 
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Figure  E.12:  Control  angle  deflections  200  km  with  yaw  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.13:  Angular  displacement  at  200km  with  roll,  pitch 
and  yaw  offset  10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.14:  Angular  rate  at  200km  with  roll,  pitch  and  yaw 
offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.15:  Angular  Acceleration  at  200km  with  roll,  pitch 
and  yaw  offset  10°. 


Control  Angle  Deflections 


Figure  E.16:  Control  angle  deflections  200  km  with  roll,  pitch 
and  yaw  offset  10°. 
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rad/s 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.17:  Angular  displacement  at  300km  with  roll  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.18:  Angular  rate  at  300km  with  roll  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.19:  Angular  Acceleration  at  300km  with  roll  offset 
10°. 
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Figure  E.20:  Control  angle  deflections  300  km  with  roll  offset 
10°. 
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rad/s 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.21:  Angular  displacement  at  300km  with  pitch  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.22:  Angular  rate  at  300km  with  pitch  offset  10°. 


112 


Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.23:  Angular  Acceleration  at  300km  with  pitch  offset 
10°. 


Control  Angle  Deflections 


Figure  E.24:  Control  angle  deflections  300  km  with  pitch  off¬ 
set  10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.25:  Angular  displacement  at  300km  with  yaw  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.26:  Angular  rate  at  300km  with  yaw  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.27:  Angular  Acceleration  at  300km  with  yaw  offset 
10°. 
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Figure  E.28:  Control  angle  deflections  at  300km  with  yaw 

offset  10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.29:  Angular  displacement  at  300km  with  roll,  pitch 
and  yaw  offset  10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.30:  Angular  rate  at  300km  with  roll,  pitch  and  yaw 
offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.31:  Angular  Acceleration  at  300km  with  roll,  pitch 
and  yaw  offset  10°. 
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Figure  E.32:  Control  angle  deflections  at  300km  with  roll, 

pitch  and  yaw  offset  10°. 
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rad/s 


Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.33:  Angular  displacement  at  400km  with  roll  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.34:  Angular  rate  at  400knr  with  roll  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.35:  Angular  Acceleration  at  400km  with  roll  offset 
10°. 
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Figure  E.36:  Control  angle  deflections  400  km  with  roll  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.37:  Angular  displacement  at  400km  with  pitch  offset 
10°. 


Time,  sec  x  1 

Figure  E.38:  Angular  rate  at  400km  with  pitch  offset  10°. 
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Figure  E.39:  Angular  Acceleration  at  400km  with  pitch  offset 
10°. 
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Figure  E.40:  Control  angle  deflections  400  km  with  pitch  off¬ 
set  10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.41:  Angular  displacement  at  400km  with  yaw  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.42:  Angular  rate  at  400km  with  yaw  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.43:  Angular  Acceleration  at  400km  with  yaw  offset 
10°. 
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Figure  E.44:  Control  angle  deflections  400  km  with  yaw  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.45:  Angular  displacement  at  400km  with  roll,  pitch 
and  yaw  offset  10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.46:  Angular  rate  at  400km  with  roll,  pitch  and  yaw 
offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.47:  Angular  Acceleration  at  400km  with  roll,  pitch 
and  yaw  offset  10°. 
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Control  Angle  Deflections 


Figure  E.48:  Control  angle  deflections  400  km  with  roll,  pitch 
and  yaw  offset  10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.49:  Angular  displacement  at  500km  with  roll  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.50:  Angular  rate  at  500km  with  roll  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.51:  Angular  Acceleration  at  500km  with  roll  offset 
10°. 
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Figure  E.52:  Control  angle  deflections  500  km  with  roll  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.53:  Angular  displacement  at  500km  with  pitch  offset 
10°. 


Time,  sec  x  1 

Figure  E.54:  Angular  rate  at  500km  with  pitch  offset  10°. 
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Figure  E.55:  Angular  Acceleration  at  500km  with  pitch  offset 
10°. 
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Figure  E.56:  Control  angle  deflections  500  km  with  pitch  off¬ 
set  10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.57:  Angular  displacement  at  500km  with  yaw  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.58:  Angular  rate  at  500km  with  yaw  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.59:  Angular  Acceleration  at  500km  with  yaw  offset 
10°. 
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Figure  E.60:  Control  angle  deflections  500  km  with  yaw  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 


Figure  E.61:  Angular  displacement  at  500km  with  roll,  pitch 
and  yaw  offset  10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.62:  Angular  rate  at  500km  with  roll,  pitch  and  yaw 
offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.63:  Angular  Acceleration  at  500km  with  roll,  pitch 
and  yaw  offset  10°. 
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Figure  E.64:  Control  angle  deflections  500  km  with  roll,  pitch 
and  yaw  offset  10°. 
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Figure  E.65:  Angular  displacement  at  600km  with  roll  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.66:  Angular  rate  at  600km  with  roll  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.67:  Angular  Acceleration  at  600km  with  roll  offset 
10°. 
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Figure  E.68:  Control  angle  deflections  600  km  with  roll  offset 
10°. 
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Figure  E.69:  Angular  displacement  at  600km  with  pitch  offset 
10°. 


Time,  sec  x1( 

Figure  E.70:  Angular  rate  at  600km  with  pitch  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 
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Figure  E.71:  Angular  Acceleration  at  600km  with  pitch  offset 
10°. 
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Figure  E.72:  Control  angle  deflections  600  km  with  pitch  off¬ 
set  10°. 
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Figure  E.73:  Angular  displacement  at  600km  with  yaw  offset 
10°. 


Roll,  Pitch,  &  Yaw  Angular  Rate  vs.  Time 


Figure  E.74:  Angular  rate  at  600km  with  yaw  offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.75:  Angular  Acceleration  at  600km  with  yaw  offset 
10°. 
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Figure  E.76:  Control  angle  deflections  600  km  with  yaw  offset 
10°. 
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Roll,  Pitch,  &  Yaw  Position  vs.  Time 
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Figure  E.77:  Angular  displacement  at  600km  with  roll,  pitch 
and  yaw  offset  10°. 
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Figure  E.78:  Angular  rate  at  600km  with  roll,  pitch  and  yaw 
offset  10°. 
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Roll,  Pitch,  &  Yaw  Angular  Acceleration  vs.  Time 


Figure  E.79:  Angular  Acceleration  at  600km  with  roll,  pitch 
and  yaw  offset  10°. 
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Figure  E.80:  Control  angle  deflections  600  km  with  roll,  pitch 
and  yaw  offset  10°. 
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Appendix  F.  Original  Design  Nonlinear  Drag  Equations 

F.l  External  Torques  from  Atmospheric  Drag  (Nonlinear) 

Satellites  in  LEO  (LEO)  (between  130  km  to  600  km)  experience  an  aerody¬ 
namic  drag  force  which  is  given  by: 

Fdrag  =  -pV^C'/^O  (6.1) 

Where  Fdrag  =drag  force  (N),  p  =atmospheric  density  (kg/rn3),  V  =velocity  (m/s), 
Cn  =  2  drag  coefficient  and  A  =projected  area  (nr2).  Altitudes  above  125  km 
altitude  are  in  the  free  molecular  flow  regime  [1,  page  318].  In  the  free  molecular 
flow  regime,  particles  are  typically  modeled  either  as  specular  or  diffuse  reflections.  A 
specular  reflection  assumes  that  molecules  are  perfectly  elastic  where  the  tangential 
velocity  is  constant  and  the  normal  velocity  is  reversed.  The  diffuse  model  assumes 
the  molecules  are  reflected  in  a  diffuse  manner  and  have  no  memory  of  previous 
velocities.  Either  model  imparts  a  force  normal  to  the  surface.  Spacecraft  in  LEO 
will  always  experience  a  small  drag  torque  since  it’s  not  possible  to  locate  the  CG  to 
the  exact  geometric  center  and  therefore  will  need  some  way  to  control  it’s  attitude. 
Assumptions  that  were  made  for  the  following  drag  torque  equations  are: 

•  The  drag  panels  are  only  visible  to  the  incoming  wind  only  on  the  front/out- 
ward  facing  sides. 

•  The  drag  coefficient  was  chosen  to  be  between  1.5  <  Cd  <  2.5  . 

•  Atmospheric  density  is  constant  and  averaged  over  the  orbit. (neglecting  solar 
effects  and  atmospheric  perturbations). 

F.1.1  Drag  Torque  Equations.  The  basic  design  being  modeled  as  a  cube 
shaped  spacecraft  that  has  arms  with  4  drag  panels  controlling  pitch,  roll  and  yaw. 
To  pitch  up  or  down,  the  top  or  bottom  panel  will  be  swung  out  into  the  air  stream. 
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For  yaw,  the  left  and  right  panels  are  used.  To  roll  left  and  right,  all  panels  will  be 
extended  and  rotated  with  respect  to  the  arms  in  a  “propeller”  configuration.  To 
determine  the  drag  effects  on  both  the  cube  and  the  panels,  the  velocity  must  be 
known  and  is  given  by  [11,  page  70] 

V  =  jkFT  (6.2) 

V  ^ orbit 

where  the  standard  gravitational  parameter  p  =  398600  and  the  orbital  radius 
rorbit  =  6378.135 km  +  Altitude ,  and  the  velocity  vector  is  in  the  —  aq  direction.  The 
density  must  also  be  determined  and  based  on  the  altitude  range  being  modeled,  the 
density  p  varies  between  10~8  to  lO”1^  for  LEO  orbits. 

F.1.2  Drag  Effects  from  Spacecraft  Body.  The  spacecraft  body  doesn’t 
produce  any  torques  unless  the  CG  is  not  at  the  geometric  center.  In  this  model,  the 
CG  location  can  be  moved  off  center  to  produce  small  torques  for  modelling  purposes. 
First  the  roll,  pitch,  yaw  rotation  matrix  is  used  to  go  from  the  A  frame  to  the  B 
frame  as  in  equation  3.2.  The  angle  of  the  incoming  molecules  and  the  inward  normal 
directions  n  for  the  cube  surfaces  can  be  found  by  using  the  dot  product  of  the  two 
vectors  and  since  we  are  dealing  with  unit  vectors,  cos(a)  =  v  ■  n.  Since  the  inward 
normal  vectors  are  aligned  with  the  principal  axes,  and  from  equation  3.6  the  dot 
products  for  each  cube  surface  become: 


cos(aLeft)  = 

(  —  Ol  '  fyj)  — 

tdBA 

-n  (2,1) 

(6.3a) 

COS  (ot  Right) 

(—a!  •  —62) 

T>BA 

K  (2,1) 

(6.3b) 

cos  (aBot)  =  ( 

-ai  '  -63)  : 

RBA 

-  n  (3,1) 

(6.3c) 

cos(aTop)  =  ( 

co 

-C> 

1 

T3ba 

-B  (3,1) 

(6.3d) 

COS  (c^Front) 

(-ai  ■  -bi) 

r?BA 

=  K  (1,1) 

(6.3e) 
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cos (aBack)  =  (-fli  •  h)  =  -RBA{  1,1) 


(6.3f) 


To  find  the  drag  force  on  each  exposed  surface,  the  projected  area  must  be  known 
and  is  found  by  Ap  =  ff  II(cos(a))  cos(a)dA  [2,  page  251].  where  H  is  the  Heaviside 
function  and  cos(a)  is  from  equations  6.3,  a-f.  The  Heaviside  function  is  used  so 
that  only  the  areas  the  see  the  incoming  molecules  are  taken  into  account.  Since 
the  areas  of  the  faces  are  the  side  of  a  cube,  the  integral  becomes  L2,  where  L  is  the 
length  of  the  side  of  the  spacecraft  body.  The  projected  areas  for  each  side  become: 


ALeft  =  L2H( cos(aLeft))  cos(aLeft)  (6.4a) 

A-Right  L  H (cos(ot  Right))  cos Right)  (6.4b) 

ABot  =  L2H(cos(aBot))  cos(aBot)  (6.4c) 

ATop  =  L2 H (cos(aTop))  cos(aTop)  (6.4d) 

ABront  L  H (cOs(a> Front))  COS [otpront)  (6.4e) 

ABack  =  L2 H (cos (a Back))  cos(aBack)  (6.4f) 


By  using  equation  6.1  and  substituting  in  equations  6.13,  the  drag  forces  for  each 
face  in  the  A  frame  become: 


F Left  A  =  ^CDALeftpV2v  (6.5a) 

F Right  a  ~^FDARightpV  V  (6.5b) 

FBotA  —  2 CBABotpV2v  (6.5c) 

Ft  oP  A  =  ^ CDAToppV2v  (6.5d) 

FFrontA  ^ F sAprontfAd  V  (6.5e) 
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F Back  A  D-A-BackpV  V 


(6.5f) 


By  using  the  rotation  matrix  3.2  the  force  in  the  B  frame  can  be  calculated 
by  substituting  equations  6.5  and  3.2  into  the  following  equations  (note,  we  are  only 
interested  in  the  inward  normal  force). 


F, 


Lefts 


0  1  0 


rbaf> 


Left  a 


F Rights 


0  1  0 


tdBA  jp 
Ft  B  Sight  A 


FBotB  — 


0  0  1 


rbaf , 


BotA 


F 


Tops 


0  0  1 


RBAFr 


Top  A 


Front B 


1  o  o 


rbaf \ 


Front  A 


F 


Back  B 


1  o  o 


tdBA  jp 
Ft  F  Back  A 


(6.6a) 


(6.6b) 


(6.6c) 


(6.6d) 


(6.6e) 


(6.6f) 


Once  the  drag  force  is  calculated  in  the  B  frame,  the  torques  can  be  calculated 
by  T  =  rxF  [2,  page  251]  where  r  is  the  vector  from  the  CG  to  the  center  of  pressure 
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for  each  surface  and  F  is  the  drag  force  in  the  B  frame.  For  the  cube  surfaces,  r 
becomes: 


0 


rLeft 


-0.5  L 
0 


CG 


(6.7a) 


T  Right 


0 

0.5L 

0 


-  CG 


(6.7b) 


rBot 


0 

0 

0.5L 


-  CG 


(6.7c) 


T Top 


0 

0 


-  CG 


(6.7d) 


-0.5L 


1  Front 


0.5  L 
0 
0 


-  CG 


(6.7e) 


T  Back 


-0.5  L 


0 

0 


-  CG 


(6.70 


Substituting  in  equations  6.6  and  6.17,  the  torque  equations  for  each  surface 
of  the  spacecraft  body  in  the  B  frame  become: 


M Left  T  Left^^Left  B 


(6.8a) 
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Right  1  Right^-F Right B 


(6.8b) 


Msot  ~  ^Bot^FsotB 

(6.8c) 

Mj1  Qp  rT  opxFT  Qp  g 

(6.8d) 

M p  ront  T  Front^F F  ront  B 

(6.8e) 

MBack  T'BackX-FBackB 

(6.8f) 

And  finally  the  total  torques  (equations  6.8  a-f)  from  the  spacecraft  body  are  added 
to  get: 

MsatBody  =  MLeft  +  MRight  +  MBot  +  MTop  +  MFront  +  MBack  (6.9) 

F.1.3  Drag  Effects  from  Spacecraft  Drag  Panels.  By  going  through  a  similar 
process,  the  torques  from  the  drag  panels  can  be  found.  The  drag  panels  adds  some 
complexity  to  the  problem  since  there  are  two  sets  of  rotations,  one  going  from  the 
A  frame  to  the  B  frame  (equation  3.2)  and  the  other  going  from  the  B  frame  to  the 
F  frame  which  by  using  the  roll,  pitch,  yaw  rotation  sequence  results  in: 


- 1 

o 

o 

c(02)  0  -s(02) 

c{4>  3)  s(03)  0 

0  c(0i)  s(4>  i) 

0  1  0 

-s(03)  c(03)  0 

0  -s(0i)  c((f)  i) 

s(02)  0  c(02) 

0  0  1 

(6.10) 


To  get  from  the  A  frame  to  the  F  frame  for  each  drag  panel,  R FA  =  RFBRBA 
where  RBA  is  from  equation  3.2  and  each  panel  has  its  specific  rotation  matrix  and 
set  of  angles  from  the  B  frame  to  the  F  frame: 
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rfb 


Top 


R 


FB 


Bot 


1 

0 

0 

0 

c{<t>TPan ) 

s(4>TPan ) 

0 

—  s(4>TPan) 

c(4>TPan) 

- 

1 

0 

0 

0 

c{(f>BPan ) 

s(4>BPan) 

0 

S^&BPan) 

C^&BPan) 

1 

0 

0 

c(  ( pTArm )  0  4*TArm) 

0  1  0 

&(  (pTArm)  0  c(  (f) T  Arm ) 


c(0) 

«(0) 

0 

-s(0) 

c(0) 

0 

0 
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c(0)  s(0)  0 

-s(0)  c(0)  0 

0  0  1 
(6.11b) 

c(0LArm)  ^(^Ly Irm)  0 

s(.4*LArm)  ^■i^LArm)  0 


0 

(.R*  LPan) 

c(4>LPan ) 

5(0) 

0  c(0) 

0 

0  1 

(6.11c) 
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0 
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c(0) 

0  — s(0) 

c(  ( pRArm ) 

s(  (jjRArm ) 

0 

tdFB  _ 

Ft  Right 

0 

c{(j)RPan ) 

s{(f>RPan ) 

0 

1  0 

S{  ( pRArm ) 

c(  < pRArm ) 

0 

0 

—  s{(f)RPan) 

c(4>RPan) 

m 

0  c(0) 

0 

0 

1 

(6. lid) 


where  (j)TPad,  (pBPad,  4>LPad ,  and  (j)RPad,  are  the  top,  bottom,  left  and  right  drag  panel 
angles  used  for  roll  control  and  range  from  ±45°.  <j)TArmi  4> BArm ,  4>LArm,  and  (pRArmi 
are  the  top,  bottom,  left  and  right  arms  that  are  connected  to  the  drag  panels.  The 
angles  between  the  b1  direction  and  the  arms  range  from  0  —  90°.  In  equations  6.11a 
and  6. lid,  the  arm  angles  have  a  negative  sign  so  that  all  angles  will  be  positive  so 
the  controller  will  only  have  to  input  positive  angles.  One  of  the  rotation  angles  in 
each  equation  are  always  zero,  this  is  due  to  the  fact  that  there  are  only  two  degrees 
of  freedom  for  each  panel. 


The  rest  of  the  torque  equations  follow  the  same  process  as  before  (see  equa¬ 
tions  6.3),  the  angles  between  the  incoming  molecules  and  the  inward  normal  direc¬ 
tions  are 
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cos  (a  Left)  =  (-ai  •  f2)  =  -RFA(2,1) 

(6.12a) 

cos  (a  Right)  =  (-ai  •  -f2)  =  RFA(2,1) 

(6.12b) 

cos (aBot)  =  (-ai  ■  ~h)  =  RFA( 3,1) 

(6.12c) 

cos(aTop)  =  (-ai  •  /3)  = 

(6.12d) 

The  projected  areas  from  the  panels  become: 

ALeft  =  L2 H (cos(aLeft))  cos (aLeft) 

(6.13a) 

-A-Right  L  H (cos(cr mghtf)  cos 

(6.13b) 

ABot  =  L2H(cos(aBot))  cos(aBot) 

(6.13c) 

ATop  =  L2 H (cos(aTop))  cos (aTop) 

(6.13d) 

where  L  is  the  side  length  of  each  square  drag  panel. 

Calculate  forces  from  the  drag  panels  in  the  A  frame: 

F Left  a  ~  ~jpDABeftpV2V 

(6.14a) 

FRightA  7 jFJ £)Afiightp\ 

(6.14b) 

FBotA  =  -jpDABotpV2v 

(6.14c) 

FTova  =  aFdAtopPF2v 

(6.14d) 
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Calculate  force  from  drag  panels  in  the  F  frame  and  pull  out  the  inward  normal 


component: 


F, 


Left  p 


0  1  0 


p>F  A  -p 

-K'  Left-T  i 


Left  "  Left  a 


F] 


Right  p 


o  1  o 


tdF  A  ZT1 

■‘A'  Right  F  Right  a 


FBotF  ~ 


0  0  1 


V)F  A  77 

FL  BotF  BotA 


Ft  op  p 


0  0  1 


P>F  A  771 
it  TYm-fg 


Topr  Top  a 


(6.15a) 


(6.15b) 


(6.15c) 


(6.15d) 


Convert  the  inward  normal  component  of  each  surface  in  the  F  frame  to  the 
B  frame. 


Fjjeftp  =  (RFBLeft)TFLeftF  (6.16a) 

FmghtB  =  (RFB  Right)T  F Right F  (6.16b) 

FBotB  —  (RFB  Bot)T  FpotF  (6.16c) 

Ftopb  =  {RFB Top)T  Fpopp  (6.16d) 


Find  r  from  the  CG  to  the  center  of  pressure  for  each  panel: 
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-  CG 


(6.17a) 


rLeft 


0.5 L  (LiArm  -f-  0.5Lpanei)  COS [(pLArm) 

0.5L  (Z/J4rm  4~  0.5 Lpane{^  sill ((pLArm) 

0 


T Right 


0.5 L  [LiArm  “t-  0.5 Lpanel)  COS^CpftArm) 
0.5 L  -f-  ( L/Arm  "F  0.5 Lpanel)  sill((/)/jy4rm) 


0 


(6.17b) 


rBot 


0.5L  ( LArm  "F  0.5 Lpanel)  COs(0 p Arm) 


0 


-CG 


0.5 L  ( L/Arm  ~f"  0.5 Lpanel)  Sll\((p  B  Arm) 


(6.17c) 


^ Top 


0.5L  ( LArm  ~b  0.5 Lpanel)  COS^CpT Arm) 


0 


-CG 


0.5L  (L/Arm  “b  0 .5Lpanei')  sill(07’y4rm) 


(6.17d) 


where  L  is  the  side  length  of  the  spacecraft  body,  Lat m  is  the  length  of  the  control 
arm,  LPanei  is  the  side  length  of  the  square  panels,  (pLArm,  (pRArm ,  (pBArm  and  (pTArm 
are  the  angles  of  the  arms  from  their  retracted  positions. 


Finally,  the  torque  equations  for  the  drag  panels  with  all  the  substitutions 
become: 


Ft  Left  T Left^-F^eft ^ 

(6.18a) 

M Right  ^  Right^-F Right  b 

(6.18b) 

M-Bot  T Bot^^BotB 

(6.18c) 

Mxop  ^Top^Ftopb 

(6.18d) 

(6.18e) 
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And  finally  the  total  torques  (equations  6.8  a-f)  from  the  spacecraft  body  are 
added  to  get: 

M Panels  =  MLeft  +  MRight  +  MBot  +  MTop  (6.19) 

And  the  total  torques  are: 


M Total  -MsatBody  A  M Panels 


(6.20) 
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